Commits

Anonymous committed 93d1df7

Some more work on sink particles and the idealized MHD collapse problems.
Merged particles which are of type MUST_REFINE are once again deleted.
This stopped a problem ... is it creating new ones? Libby?

Comments (0)

Files changed (6)

src/enzo/Grid_ReturnParticleEntry.C

     for (int j = 0; j < Size; j++) {
       if (id == List[j].Number) {
 	if (Flag[j] >= 0) {
+	  if (ParticleType[i] ==  PARTICLE_TYPE_MUST_REFINE) 
+	    ParticleMass[i] = FLOAT_UNDEFINED;
+	  else 
+	    ParticleMass[i] = tiny_number;	  
 	  ParticleType[i] = PARTICLE_TYPE_DARK_MATTER;
-	  ParticleMass[i] = tiny_number;
 	}
-//	if (Flag[j] >= 0) {
-//	  ParticleMass[i] = FLOAT_UNDEFINED;
-	break;
       }
+      
     }
-
+    
   }
-
 }

src/enzo/Grid_StarParticleHandler.C

 
     /* This creates sink particles which suck up mass off the grid. */
 
-    //if (STARMAKE_METHOD(SINK_PARTICLE))     printf("   Sink Particle\n"); 
+    //    if (STARMAKE_METHOD(SINK_PARTICLE))     printf("   Sink Particle\n"); 
     //if (level == MaximumRefinementLevel)     printf("   Max Refinement\n"); 
     if (STARMAKE_METHOD(SINK_PARTICLE) && level == MaximumRefinementLevel) {
       /* Set the density threshold by using the mass in a cell which
 	if (CellFlaggingMethod[method] == 2)
 	  SinkParticleMassThreshold = MinimumMassForRefinement[method]*
 	    pow(RefineBy, level*MinimumMassForRefinementLevelExponent[method]);
-	if (CellFlaggingMethod[method] == 6)
+	if (CellFlaggingMethod[method] == 6) { 
 	  JeansLengthRefinement = RefineByJeansLengthSafetyFactor;
+	}
       }
 
       if(BigStarFormation){

src/enzo/ReadStarParticleData.C

   if (!AddParticleAttributes) {
     if (fscanf(fptr, "NumberOfStarParticles = %"ISYM"\n",
 	       &NumberOfStarParticles) != 1) {
-      ENZO_FAIL("Error reading NumberOfStarParticles.\n");
+      //      ENZO_FAIL("Error reading NumberOfStarParticles.\n");
 
     }
   } else 

src/enzo/hydro_rk/CollapseMHD3DInitialize.C

 /
 /
 ************************************************************************/
+#ifdef USE_MPI
+#include <mpi.h>
+#endif /* USE_MPI */
 
 #include <string.h>
 #include <stdio.h>
 #include "Hierarchy.h"
 #include "LevelHierarchy.h"
 #include "TopGridData.h"
+#include "CommunicationUtilities.h"
 
 void WriteListOfFloats(FILE *fptr, int N, float floats[]);
 void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
 		      float *VelocityUnits, FLOAT Time);
 
 int CollapseMHD3DInitialize(FILE *fptr, FILE *Outfptr, 
-			    HierarchyEntry &TopGrid, TopGridData &MetaData)
+			    HierarchyEntry &TopGrid, TopGridData &MetaData, int SetBaryonFields)
 {
-  char *DensName = "Density";
-  char *TEName   = "TotalEnergy";
-  char *GEName   = "GasEnergy";
-  char *Vel1Name = "x-velocity";
-  char *Vel2Name = "y-velocity";
-  char *Vel3Name = "z-velocity";
-  char *BxName = "Bx";
-  char *ByName = "By";
-  char *BzName = "Bz";
-  char *PhiName = "Phi";
-  char *DebugName = "Debug";
-  char *Phi_pName = "Phip";
-  char *GravPotenName = "PotentialField";
-  char *Acce1Name = "AccelerationField1";
-  char *Acce2Name = "AccelerationField2";
-  char *Acce3Name = "AccelerationField3";
+  const char *DensName = "Density";
+  const char *TEName   = "TotalEnergy";
+  const char *GEName   = "GasEnergy";
+  const char *Vel1Name = "x-velocity";
+  const char *Vel2Name = "y-velocity";
+  const char *Vel3Name = "z-velocity";
+  const char *ElectronName = "Electron_Density";
+  const char *HIName    = "HI_Density";
+  const char *HIIName   = "HII_Density";
+  const char *HeIName   = "HeI_Density";
+  const char *HeIIName  = "HeII_Density";
+  const char *HeIIIName = "HeIII_Density";
+  const char *HMName    = "HM_Density";
+  const char *H2IName   = "H2I_Density";
+  const char *H2IIName  = "H2II_Density";
+  const char *DIName    = "DI_Density";
+  const char *DIIName   = "DII_Density";
+  const char *HDIName   = "HDI_Density";
+  const char *BxName = "Bx";
+  const char *ByName = "By";
+  const char *BzName = "Bz";
+  const char *PhiName = "Phi";
+  const char *DebugName = "Debug";
+  const char *Phi_pName = "Phip";
+  const char *GravPotenName = "PotentialField";
+  const char *Acce1Name = "AccelerationField1";
+  const char *Acce2Name = "AccelerationField2";
+  const char *Acce3Name = "AccelerationField3";
 
   /* declarations */
 
   int UseParticles    = FALSE;
   float MediumDensity = 1.0, 
     MediumPressure = 1.0;
+  float MediumTemperature = FLOAT_UNDEFINED;
   int   SphereType[MAX_SPHERES];
   float SphereDensity[MAX_SPHERES],
     SpherePressure[MAX_SPHERES],
     SphereCutOff[MAX_SPHERES],
     SphereAng1[MAX_SPHERES],
     SphereAng2[MAX_SPHERES];
+
   int	SphereNumShells[MAX_SPHERES];
   FLOAT SphereRadius[MAX_SPHERES],
     SphereCoreRadius[MAX_SPHERES],
     SpherePosition[MAX_SPHERES][MAX_DIMENSION];
   float Bnaught = 0.0;
   float theta_B = 0.0;
+  int Bdirection = 2;
 
   for (sphere = 0; sphere < MAX_SPHERES; sphere++) {
     SphereRadius[sphere]     = 1.0;
     ret += sscanf(line, "UseParticles = %"ISYM, &UseParticles);
     ret += sscanf(line, "MediumDensity = %"FSYM, &MediumDensity);
     ret += sscanf(line, "MediumPressure = %"FSYM, &MediumPressure);
+    ret += sscanf(line, "MediumTemperature = %"FSYM, &MediumTemperature);
     ret += sscanf(line, "UniformVelocity = %"FSYM" %"FSYM" %"FSYM, 
 		  UniformVelocity, UniformVelocity+1,
 		  UniformVelocity+2);
     ret += sscanf(line, "InitialBField = %"FSYM, &Bnaught);
     ret += sscanf(line, "theta_B = %"FSYM, &theta_B);
- 
+    ret += sscanf(line, "Bdirection = %"ISYM, &Bdirection);
+
     if (sscanf(line, "SphereType[%"ISYM"]", &sphere) > 0)
       ret += sscanf(line, "SphereType[%"ISYM"] = %"ISYM, &sphere,
 		    &SphereType[sphere]);
     double rhoc = ksi_e*ksi_e*f*cs*cs/(re*re*4*M_PI*G);
 
     SphereDensity[0] = rhoc;
-    //    MediumDensity = rhoc/14.0;
+    MediumDensity *= rhoc/14.0; // in this case medium density is the ratio to decrease the outside density 
     
-    MediumPressure = rhoc*cs*cs/14.0;
+    MediumPressure = rhoc/14.0 *cs*cs/Gamma;
     double m_be = pow(f,1.5)*1.18*pow(cs,4)/pow(G,1.5)/sqrt(MediumPressure);
     double msun = 1.989e33;
     m_be /= msun;
 
 
   MediumDensity /= rhou;
-  MediumPressure /= presu;
+  if (MediumTemperature > 0) {
+    float MediumEnergy;
+    float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, 
+      VelocityUnits;
+    float tmp1, tmp2, tmp3;
+    GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, 
+	     &VelocityUnits, MetaData.Time);
+    MediumEnergy = MediumTemperature / TemperatureUnits / ((Gamma-1.0)*Mu);
+    MediumPressure = (Gamma-1.0) * MediumDensity * MediumEnergy;
+  } else {
+    MediumPressure /= presu;
+  }
 
   Bnaught /= bfieldu;
 
 
   printf("rhoc=%"GSYM", rhom=%"GSYM", pm=%"GSYM"\n", SphereDensity[0], MediumDensity, MediumPressure);
 
-
-  if (TopGrid.GridData->CollapseMHD3DInitializeGrid(
-	     n_sphere, SphereRadius,
-	     SphereCoreRadius, SphereDensity,
-	     SpherePressure, SphereSoundVelocity, SpherePosition, 
-	     SphereAngVel, Bnaught, theta_B, 
-	     SphereType,
-             MediumDensity, MediumPressure, 0) == FAIL) {
-    fprintf(stderr, "Error in CollapseTestInitializeGrid.\n");
-    return FAIL;
+  HierarchyEntry *CurrentGrid; // all level 0 grids on this processor first
+  CurrentGrid = &TopGrid;
+  int count = 0;
+  while (CurrentGrid != NULL) {
+    printf("count %i %i\n", count++, MyProcessorNumber);
+    if (CurrentGrid->GridData->CollapseMHD3DInitializeGrid(
+						      n_sphere, SphereRadius,
+						      SphereCoreRadius, SphereDensity,
+						      SpherePressure, SphereSoundVelocity, SpherePosition, 
+						      SphereAngVel, SphereTurbulence, Bnaught, theta_B, Bdirection,
+						      SphereType,
+						      MediumDensity, MediumPressure, 0, SetBaryonFields) == FAIL) {
+      fprintf(stderr, "Error in CollapseMHD3DInitializeGrid.\n");
+      return FAIL;
+    }
+    CurrentGrid = CurrentGrid->NextGridThisLevel;
   }
 
-  /* Convert minimum initial overdensity for refinement to mass
-     (unless MinimumMass itself was actually set). */
 
-  if (MinimumMassForRefinement[0] == FLOAT_UNDEFINED) {
-    MinimumMassForRefinement[0] = MinimumOverDensityForRefinement[0];
-    for (int dim = 0; dim < MetaData.TopGridRank; dim++)
-      MinimumMassForRefinement[0] *=(DomainRightEdge[dim]-DomainLeftEdge[dim])/
-	float(MetaData.TopGridDims[dim]);
-  }
+  if (SetBaryonFields) {
 
-  /* If requested, refine the grid to the desired level. */
-
-  if (RefineAtStart) {
-
-    /* Declare, initialize and fill out the LevelArray. */
-
-    LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY];
-    for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
-      LevelArray[level] = NULL;
-    AddLevel(LevelArray, &TopGrid, 0);
-
-    /* Add levels to the maximum depth or until no new levels are created,
-       and re-initialize the level after it is created. */
-
-    for (level = 0; level < MaximumRefinementLevel; level++) {
-
-      if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) {
-	fprintf(stderr, "Error in RebuildHierarchy.\n");
-	return FAIL;
-      }
-      if (LevelArray[level+1] == NULL)
-	break;
-
-      LevelHierarchyEntry *Temp = LevelArray[level+1];
-      while (Temp != NULL) {
-	if (Temp->GridData->CollapseMHD3DInitializeGrid(
-				n_sphere, SphereRadius,
-				SphereCoreRadius, SphereDensity,
-				SpherePressure, SphereSoundVelocity, SpherePosition, SphereAngVel, Bnaught,theta_B,
-				SphereType, MediumDensity, MediumPressure, level+1) == FAIL) {
-	  fprintf(stderr, "Error in Collapse3DInitializeGrid.\n");
+    // Compute Velocity Normalization
+    double v_rms  = 0;
+    double Volume = 0;
+    Eflt fac = 1;
+    
+    if (SphereTurbulence[0] > 0.) {
+      CurrentGrid = &TopGrid;
+      while (CurrentGrid != NULL) {
+	if (CurrentGrid->GridData->PrepareVelocityNormalization(&v_rms, &Volume) == FAIL) {
+	  fprintf(stderr, "Error in PrepareVelocityNormalization.\n");
 	  return FAIL;
 	}
-	Temp = Temp->NextGridThisLevel;
+	CurrentGrid = CurrentGrid->NextGridThisLevel;
+	fprintf(stderr, "Prepared: v_rms, Volume: %"GSYM"  %"GSYM"\n", v_rms, Volume);
       }
-    } // end: loop over levels
-
-    /* Loop back from the bottom, restoring the consistency among levels. */
-
-    for (level = MaximumRefinementLevel; level > 0; level--) {
-      LevelHierarchyEntry *Temp = LevelArray[level];
-      while (Temp != NULL) {
-	if (Temp->GridData->ProjectSolutionToParentGrid(
-				   *LevelArray[level-1]->GridData) == FAIL) {
-	  fprintf(stderr, "Error in grid->ProjectSolutionToParentGrid.\n");
+      
+#ifdef USE_MPI
+      CommunicationAllReduceValues(&v_rms, 1, MPI_SUM);
+      CommunicationAllReduceValues(&Volume, 1, MPI_SUM);
+#endif
+      fprintf(stderr, "v_rms, Volume: %"GSYM"  %"GSYM"\n", v_rms, Volume);
+      // Carry out the Normalization
+      v_rms = sqrt(v_rms/Volume); // actuall v_rms
+      fac = SphereSoundVelocity[0]*SphereTurbulence[0]/v_rms;
+      
+      CurrentGrid = &TopGrid;
+      while (CurrentGrid != NULL) {
+	if (CurrentGrid->GridData->NormalizeVelocities(fac) == FAIL) {
+	  fprintf(stderr, "Error in grid::NormalizeVelocities.\n");
 	  return FAIL;
 	}
-	Temp = Temp->NextGridThisLevel;
+	CurrentGrid = CurrentGrid->NextGridThisLevel;
       }
     }
+    if (fac != 0. ) SphereTurbulence[0] = fac;
 
-  } // end: if (RefineAtStart)
+    /* Convert minimum initial overdensity for refinement to mass
+       (unless MinimumMass itself was actually set). */
 
-  /* set up field names and units */
+    if (MinimumMassForRefinement[0] == FLOAT_UNDEFINED) {
+      MinimumMassForRefinement[0] = MinimumOverDensityForRefinement[0];
+      for (int dim = 0; dim < MetaData.TopGridRank; dim++)
+	MinimumMassForRefinement[0] *=(DomainRightEdge[dim]-DomainLeftEdge[dim])/
+	  float(MetaData.TopGridDims[dim]);
+    }
 
-  int count = 0;
-  DataLabel[count++] = DensName;
-  DataLabel[count++] = Vel1Name;
-  DataLabel[count++] = Vel2Name;
-  DataLabel[count++] = Vel3Name;
-  DataLabel[count++] = TEName;
-  if (DualEnergyFormalism) {
-    DataLabel[count++] = GEName;
-  }
-  if (HydroMethod == MHD_RK) {
-    DataLabel[count++] = BxName;
-    DataLabel[count++] = ByName;
-    DataLabel[count++] = BzName;
-    DataLabel[count++] = PhiName;
-  }
+    /* If requested, refine the grid to the desired level. */
 
-  if(UseDivergenceCleaning){
-    DataLabel[count++] = Phi_pName;
-    DataLabel[count++] = DebugName;
-  }
+    if (RefineAtStart) {
 
-  if (WritePotential) {
-    DataLabel[count++] = GravPotenName;
-    DataLabel[count++] = Acce1Name;
-    DataLabel[count++] = Acce2Name;
-    DataLabel[count++] = Acce3Name;
-  }
+      /* Declare, initialize and fill out the LevelArray. */
 
-  for (i = 0; i < count; i++) {
-    DataUnits[i] = NULL;
-  }
+      LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY];
+      for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
+	LevelArray[level] = NULL;
+      AddLevel(LevelArray, &TopGrid, 0);
+
+      /* Add levels to the maximum depth or until no new levels are created,
+	 and re-initialize the level after it is created. */
+
+      for (level = 0; level < MaximumRefinementLevel; level++) {
+
+	if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) {
+	  fprintf(stderr, "Error in RebuildHierarchy.\n");
+	  return FAIL;
+	}
+	if (LevelArray[level+1] == NULL)
+	  break;
+
+	LevelHierarchyEntry *Temp = LevelArray[level+1];
+	while (Temp != NULL) {
+	  if (Temp->GridData->CollapseMHD3DInitializeGrid(
+							  n_sphere, SphereRadius,
+							  SphereCoreRadius, SphereDensity,
+							  SpherePressure, SphereSoundVelocity, 
+							  SpherePosition, SphereAngVel, SphereTurbulence,  
+							  Bnaught, theta_B, Bdirection,
+							  SphereType, MediumDensity, MediumPressure, level+1,
+							  SetBaryonFields) == FAIL) {
+	    fprintf(stderr, "Error in Collapse3DInitializeGrid.\n");
+	    return FAIL;
+	  }
+	  Temp = Temp->NextGridThisLevel;
+	}
+      } // end: loop over levels
+    
+      /* Loop back from the bottom, restoring the consistency among levels. */
+
+      for (level = MaximumRefinementLevel; level > 0; level--) {
+	LevelHierarchyEntry *Temp = LevelArray[level];
+	while (Temp != NULL) {
+	  if (Temp->GridData->ProjectSolutionToParentGrid(
+							  *LevelArray[level-1]->GridData) == FAIL) {
+	    fprintf(stderr, "Error in grid->ProjectSolutionToParentGrid.\n");
+	    return FAIL;
+	  }
+	  Temp = Temp->NextGridThisLevel;
+	}
+      }
+
+    } // end: if (RefineAtStart)
+
+    /* set up field names and units */
+
+    int count = 0;
+    DataLabel[count++] = (char*) DensName;
+    DataLabel[count++] = (char*) Vel1Name;
+    DataLabel[count++] = (char*) Vel2Name;
+    DataLabel[count++] = (char*) Vel3Name;
+    DataLabel[count++] = (char*) TEName;
+    if (DualEnergyFormalism) {
+      DataLabel[count++] = (char*) GEName;
+    }
+    if (HydroMethod == MHD_RK) {
+      DataLabel[count++] = (char*) BxName;
+      DataLabel[count++] = (char*) ByName;
+      DataLabel[count++] = (char*) BzName;
+      DataLabel[count++] = (char*) PhiName;
+    }
+
+    if (MultiSpecies) {
+      DataLabel[count++] = (char*) ElectronName;
+      DataLabel[count++] = (char*) HIName;
+      DataLabel[count++] = (char*) HIIName;
+      DataLabel[count++] = (char*) HeIName;
+      DataLabel[count++] = (char*) HeIIName;
+      DataLabel[count++] = (char*) HeIIIName;
+      if (MultiSpecies > 1) {
+	DataLabel[count++] = (char*) HMName;
+	DataLabel[count++] = (char*) H2IName;
+	DataLabel[count++] = (char*) H2IIName;
+      }
+      if (MultiSpecies > 2) {
+	DataLabel[count++] = (char*) DIName;
+	DataLabel[count++] = (char*) DIIName;
+	DataLabel[count++] = (char*) HDIName;
+      }
+    }  // if Multispecies
+
+    if(UseDivergenceCleaning){
+      DataLabel[count++] = (char*) Phi_pName;
+      DataLabel[count++] = (char*) DebugName;
+    }
+
+    if (WritePotential) {
+      DataLabel[count++] = (char*) GravPotenName;
+      DataLabel[count++] = (char*) Acce1Name;
+      DataLabel[count++] = (char*) Acce2Name;
+      DataLabel[count++] = (char*) Acce3Name;
+    }
+
+    for (i = 0; i < count; i++) {
+      DataUnits[i] = NULL;
+    }
+
+  } // end if  SetBaryonField
 
   return SUCCESS;
 

src/enzo/hydro_rk/Grid_CollapseMHD3DInitializeGrid.C

 /
 /  written by: Peng Wang
 /  date:       June, 2007
-/  modified1:
+/  modified1: Tom Abel 2010, add turbulence generator to this problem as well
 /
 /
 ************************************************************************/
 	     float *TemperatureUnits, float *TimeUnits,
 	     float *VelocityUnits, FLOAT Time);
 double Gaussian(double cs);
+void Turbulence_Generator(float **vel, int dim0, int dim1, int dim2, int ind, 
+			  float kmin, float kmax, float dk,
+			  FLOAT **LeftEdge, FLOAT **CellWidth, int seed);
+
 
 int grid::CollapseMHD3DInitializeGrid(int n_sphere,
 				      FLOAT r_sphere[MAX_SPHERES],
 				      float p_sphere[MAX_SPHERES],
 				      float cs_sphere[MAX_SPHERES],
 				      FLOAT sphere_position[MAX_SPHERES][MAX_DIMENSION],
-				      float omega_sphere[MAX_SPHERES], float Bnaught, float theta_B,
+				      float omega_sphere[MAX_SPHERES], 
+				      float turb_sphere[MAX_SPHERES], 
+				      float Bnaught, float theta_B,
+				      int Bdirection,
 				      int   sphere_type[MAX_SPHERES],
-				      float rho_medium, float p_medium, int level)
+				      float rho_medium, float p_medium, int level,
+				      int SetBaryonFields)
 {
+
+  const float HIIFraction = 1.2e-5;
+  const float HeIIFraction = 1.0e-14;
+  const float HeIIIFraction = 1.0e-17;
+  const float HMFraction = 2.0e-9;
+  const float H2IFraction = 2.0e-4;
+  const float H2IIFraction = 3.0e-14;
+  int TurbulenceSeed = 191105;
+  float *TurbulenceVelocity[3];
+  float HIFraction, HeIFraction, eFraction;
+  HIFraction = CoolData.HydrogenFractionByMass - HIIFraction;
+  if (MultiSpecies > 1)
+    HIFraction -= HMFraction + H2IFraction + H2IIFraction;
+  HeIFraction = 1.0 - CoolData.HydrogenFractionByMass - 
+    HeIIFraction - HeIIIFraction;
+  eFraction = HIIFraction + 0.25*HeIIFraction + 0.5*HeIIIFraction;
+  if (MultiSpecies > 1)
+    eFraction += 0.5*H2IIFraction - HMFraction;
+
+
   /* declarations */
 
-  int dim, i, j, k, m, field, sphere, size;
+  int dim, i, j, k, m, field, sphere, size, activesize;
 
-    int phip_num;
+  int phip_num;
   NumberOfBaryonFields = 0;
   FieldType[NumberOfBaryonFields++] = Density;
   FieldType[NumberOfBaryonFields++] = Velocity1;
     FieldType[NumberOfBaryonFields++] = DebugField;  
   }
 
+  int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, 
+    H2INum, H2IINum, DINum, DIINum, HDINum;
+
+  if (MultiSpecies) {
+    FieldType[DeNum    = NumberOfBaryonFields++] = ElectronDensity;
+    FieldType[HINum    = NumberOfBaryonFields++] = HIDensity;
+    FieldType[HIINum   = NumberOfBaryonFields++] = HIIDensity;
+    FieldType[HeINum   = NumberOfBaryonFields++] = HeIDensity;
+    FieldType[HeIINum  = NumberOfBaryonFields++] = HeIIDensity;
+    FieldType[HeIIINum = NumberOfBaryonFields++] = HeIIIDensity;
+    if (MultiSpecies > 1) {
+      FieldType[HMNum    = NumberOfBaryonFields++] = HMDensity;
+      FieldType[H2INum   = NumberOfBaryonFields++] = H2IDensity;
+      FieldType[H2IINum  = NumberOfBaryonFields++] = H2IIDensity;
+    }
+    if (MultiSpecies > 2) {
+      FieldType[DINum   = NumberOfBaryonFields++] = DIDensity;
+      FieldType[DIINum  = NumberOfBaryonFields++] = DIIDensity;
+      FieldType[HDINum  = NumberOfBaryonFields++] = HDIDensity;
+    }
+  }
+
   if (WritePotential) {
     FieldType[NumberOfBaryonFields++] = GravPotential;
     FieldType[NumberOfBaryonFields++] = AccelerationField1;
     return SUCCESS;
   }
 
+  //  for (i=0; i< NumberOfBaryonFields; i++)    BaryonField[i] = NULL;
+
+
+  printf("Grid_CollapseMHD3DInitialize: Setting up grid variables.\n");
+
   /* Units and parameters */
 
   float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1.0, 
   for (dim = 0; dim < GridRank; dim++) {
     size *= GridDimension[dim];
   }
+  activesize = 1;
+  for (dim = 0; dim < GridRank; dim++) {
+    activesize *= (GridDimension[dim] - 2*DEFAULT_GHOST_ZONES);
+  }
+
+  for (dim = 0; dim < GridRank; dim++) {
+    TurbulenceVelocity[dim] = new float[activesize];
+    for (int n = 0; n < activesize; n++) {
+      TurbulenceVelocity[dim][n] = 0.0;
+    }
+  }
 
   int count=0;
   for (field = 0; field < NumberOfBaryonFields; field++) {
   printf("rho_sphere=%"GSYM", cs_sphere=%"GSYM", rho_medium=%"GSYM", p_medium=%"GSYM"\n",
 	 rho_sphere[0], cs_sphere[0], rho_medium, p_medium);
   printf("r_sphere: %"GSYM"\n", r_sphere[0]);
+  printf("turb_sphere: %"GSYM"\n", turb_sphere[0]);
+
   // if use BE sphere, read in the BE sphere density profile
 
   char *filename = "be.dat";
 	for (dim = 0; dim < 3; dim++) {
 	  vel[dim] = 0.0;
 	}
+
 	Bx = 0.0;
 	By = 0.0;
-	Bz = Bnaught;
+	Bz = 0.0;
+
+	switch (Bdirection) {
+	case 0:
+	  Bx = Bnaught;
+	  break;
+	case 1:
+	  By = Bnaught;
+	  break;
+	case 2:
+	  Bz = Bnaught;
+	  break;
+	default:
+	  ENZO_FAIL("Bdirection must be 0,1,2");
+	}
 
 	/* Loop over spheres. */
 	for (sphere = 0; sphere < n_sphere; sphere++) {
 		n = n_bin -1;
 	      }
 	      rho = rho_sphere[sphere]*rho_be[n];
-	      eint = pow(cs_sphere[sphere], 2)/(Gamma-1.0);
+	      float p, cs, h, dpdrho, dpde;
+	      p = rho*pow(cs_sphere[sphere], 2)/Gamma;
+	      EOS(p, rho, eint, h, cs, dpdrho, dpde, EOSType, 1); 
+	      //	      eint = pow(cs_sphere[sphere], 2)/(Gamma-1.0);
 	      FLOAT omega = omega_sphere[sphere];
 	      vel[0] = -omega*ypos;
 	      vel[1] = omega*xpos;
     }
   }
 
+
+
+  /* Initialize turbulent velocity field */
+
+  if (SetBaryonFields && (turb_sphere[0] > 0.)) {
+
+    float k1, k2, dk;
+      k1 = 5;
+      k2 = 8.0;
+      dk = 1.0;
+
+    printf("Begin generating turbulent velocity spectrum...\n");
+    Turbulence_Generator(TurbulenceVelocity, 
+			 GridDimension[0]-2*DEFAULT_GHOST_ZONES, 
+			 GridDimension[1]-2*DEFAULT_GHOST_ZONES,
+			 GridDimension[2]-2*DEFAULT_GHOST_ZONES,
+			 4.0, k1, k2, dk,
+			 CellLeftEdge, CellWidth, TurbulenceSeed);    
+
+    printf("Turbulent spectrum generated\n");
+
+    float VelocityNormalization = 1;
+// for level > 0 grids the CloudMachNumber passed in is actuall the Velocity normalization factor
+    if (level > 0) VelocityNormalization = turb_sphere[0];
+    printf("Cloud Mach Number = %"GSYM" \n",turb_sphere[0]);
+    for (i = 0; i < 3; i++) {
+      for (n = 0; n < activesize; n++) {
+	TurbulenceVelocity[i][n] *= VelocityNormalization;
+      }
+    }
+
+
+    /* Set turbulent velocity field */
+  FLOAT x,y,z;
+  int igrid;
+    n = 0;
+    /*    for (k = 0; k < GridDimension[2]; k++) {
+      for (j = 0; j < GridDimension[1]; j++) {
+      for (i = 0; i < GridDimension[0]; i++, n++) { */
+	  for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) { 
+	    for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+	      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, n++) { 
+	  igrid = i + GridDimension[0]*(j+k*GridDimension[1]);
+	  x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
+	  y = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
+	  z = CellLeftEdge[2][k] + 0.5*CellWidth[2][k];
+	  
+	  BaryonField[ivx][igrid] += TurbulenceVelocity[0][n];
+	  BaryonField[ivy][igrid] += TurbulenceVelocity[1][n];
+	  BaryonField[ivz][igrid] += TurbulenceVelocity[2][n];
+	  BaryonField[ietot][igrid] += 
+	    0.5 * (pow(TurbulenceVelocity[0][n],2) + 
+		   pow(TurbulenceVelocity[1][n],2) + 
+		   pow(TurbulenceVelocity[2][n],2));
+
+	} 
+      }
+    }    
+
+  }
+
+
+
+
+
+
   int TestBinary = 0;
   if (TestBinary == 1 && level == 0) {
 
     ParticleVelocity[0][0] = 0.0;
     ParticleVelocity[1][0] = 0.0;
     ParticleVelocity[2][0] = 0.0;
-    ParticleAttribute[0][0] = 0.01; // creation time    
-    ParticleAttribute[1][0] = 0; // dynamical time                                                                
+    ParticleAttribute[0][0] = 0.01; // creation time
+    ParticleAttribute[1][0] = 0; // dynamical time
     ParticleAttribute[2][0] = mass_p;
   }
 
+  /* Set uniform species fractions */
 
+  int index;
+  if (MultiSpecies > 0) {
+    for (k = 0, index = 0; k < GridDimension[2]; k++)
+      for (j = 0; j < GridDimension[1]; j++)
+	for (i = 0; i < GridDimension[0]; i++, index++) {
+	  BaryonField[DeNum][index] = eFraction * BaryonField[0][index];
+	  BaryonField[HINum][index] = HIFraction * BaryonField[0][index];
+	  BaryonField[HIINum][index] = HIIFraction * BaryonField[0][index];
+	  BaryonField[HeINum][index] = HeIFraction * BaryonField[0][index];
+	  BaryonField[HeIINum][index] = HeIIFraction * BaryonField[0][index];
+	  BaryonField[HeIIINum][index] = HeIIIFraction * BaryonField[0][index];
+	}
+  }
+
+  if (MultiSpecies > 1) {
+    for (k = 0, index = 0; k < GridDimension[2]; k++)
+      for (j = 0; j < GridDimension[1]; j++)
+	for (i = 0; i < GridDimension[0]; i++, index++) {
+	  BaryonField[HMNum][index] = HMFraction * BaryonField[0][index];
+	  BaryonField[H2INum][index] = H2IFraction * BaryonField[0][index];
+	  BaryonField[H2IINum][index] = H2IIFraction * BaryonField[0][index];
+	}
+  }
+
+  if (MultiSpecies > 2) {
+    for (k = 0, index = 0; k < GridDimension[2]; k++)
+      for (j = 0; j < GridDimension[1]; j++)
+	for (i = 0; i < GridDimension[0]; i++, index++) {
+	  BaryonField[DINum][index] = CoolData.DeuteriumToHydrogenRatio * 
+	    BaryonField[HINum][index];
+	  BaryonField[DIINum][index] = CoolData.DeuteriumToHydrogenRatio * 
+	    BaryonField[HIINum][index];
+	  BaryonField[HDINum][index] = CoolData.DeuteriumToHydrogenRatio * 
+	    BaryonField[H2INum][index];
+	}
+  }
+
+  printf("Grid_CollapseMHD3DInitiailize: done with this grid\n");
   return SUCCESS;
 }

src/enzo/star_maker8.C

     ny_cell[MAX_SUPERCELL_NUMBER], nz_cell[MAX_SUPERCELL_NUMBER];
 
 
-  /*printf("Star Maker 8 running - SinkMergeDistance = %g\n", SinkMergeDistance);
+  printf("Star Maker 8 running - SinkMergeDistance = %g\n", SinkMergeDistance);
   printf("Star Maker 8: massthresh=%g, jlrefine=%g\n", *massthresh,*jlrefine);
-  printf("Star Maker 8: time = %g\n", *t); */
+  printf("Star Maker 8: time = %g\n", *t); 
 
 
   /* Compute Units. */