Commits

Michael Kuhlen  committed c48f290

Moving not yet converted test problem initializations to not_yet_converted/, this time with 'hg mv'.

  • Participants
  • Parent commits 9433149

Comments (0)

Files changed (124)

File src/enzo/problem_types/AdiabaticExpansionInitialize.C

-/***********************************************************************
-/
-/  INITIALIZE AN ADIABATIC EXPANSION TEST
-/
-/  written by: Greg Bryan
-/  date:       April, 1995
-/  modified1:
-/
-/  PURPOSE:
-/
-/  RETURNS: SUCCESS or FAIL
-/
-************************************************************************/
- 
-// This routine intializes a new simulation based on the parameter file.
-//
- 
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-#include "CosmologyParameters.h"
- 
-/* Set the mean molecular mass as in Grid_ComputeTemperatureField.C */
- 
-#define DEFAULT_MU 0.6
- 
-/* function prototypes */
- 
-void WriteListOfFloats(FILE *fptr, int N, float floats[]);
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, FLOAT Time);
- 
-int CosmologyComputeTimeFromRedshift(FLOAT Redshift, FLOAT *TimeCodeUnits);
-
-int AdiabaticExpansionInitialize(FILE *fptr, FILE *Outfptr,
-			       HierarchyEntry &TopGrid)
-{
-  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 *GPotName  = "Grav_Potential";
-
-  /* declarations */
- 
-  char line[MAX_LINE_LENGTH];
-  int ret;
- 
-  /* Error check. */
- 
-  if (!ComovingCoordinates) {
-        ENZO_FAIL("ComovingCoordinates must be TRUE!");
-  }
- 
-  /* set default parameters */
- 
-  float AdiabaticExpansionOmegaBaryonNow     = 1.0;  // standard
-  float AdiabaticExpansionOmegaCDMNow        = 0.0;  // no dark matter
-  float AdiabaticExpansionInitialTemperature = 200;  // degrees K
-  float AdiabaticExpansionInitialUniformBField[MAX_DIMENSION];  // in Gauss
-  float AdiabaticExpansionInitialVelocity    = 100;  // km/s
- 
-  float InitialVels[MAX_DIMENSION];  // Initialize arrays
-  for (int dim = 0; dim < MAX_DIMENSION; dim++) {
-    AdiabaticExpansionInitialUniformBField[dim] = 0.0;
-    InitialVels[dim] = 0.0;
-  }
-
-  /* read input from file */
- 
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
- 
-    ret = 0;
- 
-    /* read parameters */
- 
-    ret += sscanf(line, "AdiabaticExpansionOmegaBaryonNow = %"FSYM,
-		  &AdiabaticExpansionOmegaBaryonNow);
-    ret += sscanf(line, "AdiabaticExpansionOmegaCDMNow = %"FSYM,
-		  &AdiabaticExpansionOmegaCDMNow);
-    ret += sscanf(line, "AdiabaticExpansionInitialTemperature = %"FSYM,
-		  &AdiabaticExpansionInitialTemperature);
-    ret += sscanf(line, "AdiabaticExpansionInitialUniformBField = %"FSYM" %"FSYM" %"FSYM,
-		  AdiabaticExpansionInitialUniformBField,
-		  AdiabaticExpansionInitialUniformBField+1,
-		  AdiabaticExpansionInitialUniformBField+2);
-    ret += sscanf(line, "AdiabaticExpansionInitialVelocity = %"FSYM,
-		  &AdiabaticExpansionInitialVelocity);
- 
-    /* if the line is suspicious, issue a warning */
- 
-    if (ret == 0 && strstr(line, "=") && strstr(line, "AdiabaticExpansion"))
-      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
- 
-  }
-
-  /* error checking */
-  if (Mu != DEFAULT_MU) {
-    if (MyProcessorNumber == ROOT_PROCESSOR)
-      fprintf(stderr, "warning: mu =%f assumed in initialization; setting Mu = %f for consistency.\n", DEFAULT_MU);
-    Mu = DEFAULT_MU;
-  }
-
- 
-  /* Get the units so we can convert temperature later. */
- 
-  float DensityUnits=1, LengthUnits=1, TemperatureUnits=1, TimeUnits=1,
-    VelocityUnits=1, PressureUnits=1.,MagneticUnits=1., a=1,dadt=0;
-
-  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	       &TimeUnits, &VelocityUnits, InitialTimeInCodeUnits) == FAIL) {
-        ENZO_FAIL("Error in GetUnits.");
-  }
-  PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
-  MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
-
-  for (int dim = 0; dim < MAX_DIMENSION; dim++) 
-    AdiabaticExpansionInitialUniformBField[dim] /= MagneticUnits;
-  
-  /* Put inputs in a form that will be understood by InitializeUniformGrid. */
- 
-  float InitialTotalEnergy, InitialGasEnergy;
-  InitialGasEnergy = AdiabaticExpansionInitialTemperature/TemperatureUnits /
-    (Gamma - 1.0);
-  if (MultiSpecies == FALSE) InitialGasEnergy /=  DEFAULT_MU;
-  InitialTotalEnergy = InitialGasEnergy;
-  InitialVels[0] = AdiabaticExpansionInitialVelocity/VelocityUnits*1.0e5;
-  InitialTotalEnergy += 0.5*POW(InitialVels[0],2);
-  for (int dim = 1; dim < MAX_DIMENSION; dim++)
-    InitialVels[dim] = 0.0;
- 
-  /* set up grid */
- 
-  if (TopGrid.GridData->InitializeUniformGrid(
-					      AdiabaticExpansionOmegaBaryonNow,
-					      InitialTotalEnergy,
-					      InitialGasEnergy, InitialVels,
-					      AdiabaticExpansionInitialUniformBField
-					      ) == FAIL) {
-        ENZO_FAIL("Error in InitializeUniformGrid.");
-  }
- 
-  /* set up field names and units */
- 
-  int i = 0;
-  DataLabel[i++] = DensName;
-  DataLabel[i++] = TEName;
-  if (DualEnergyFormalism)
-    DataLabel[i++] = GEName;
-  DataLabel[i++] = Vel1Name;
-  DataLabel[i++] = Vel2Name;
-  DataLabel[i++] = Vel3Name;
-  if (HydroMethod == MHD_RK) {
-    DataLabel[i++] = BxName;
-    DataLabel[i++] = ByName;
-    DataLabel[i++] = BzName;
-    DataLabel[i++] = PhiName;
-  }
-  if(UseDivergenceCleaning){
-    DataLabel[i++] = Phi_pName;
-    DataLabel[i++] = DebugName;
-  }
-  if (WritePotential)
-    DataLabel[i++] = GPotName;  
- 
-
-  DataUnits[0] = NULL;
-  DataUnits[1] = NULL;
-  DataUnits[2] = NULL;
-  DataUnits[3] = NULL;
-  DataUnits[4] = NULL;
- 
-  /* Write parameters to parameter output file */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    fprintf(Outfptr, "AdiabaticExpansionOmegaBaryonNow     = %"FSYM"\n",
-	    AdiabaticExpansionOmegaBaryonNow);
-    fprintf(Outfptr, "AdiabaticExpansionOmegaCDMNow        = %"FSYM"\n",
-	    AdiabaticExpansionOmegaCDMNow);
-    fprintf(Outfptr, "AdiabaticExpansionInitialTemperature = %"FSYM"\n",
-	    AdiabaticExpansionInitialTemperature);
-    fprintf(Outfptr, "AdiabaticExpansionInitialUniformBField = ");
-    WriteListOfFloats(Outfptr, 3, AdiabaticExpansionInitialUniformBField);
-
-    fprintf(Outfptr, "AdiabaticExpansionInitialVelocity    = %"FSYM"\n\n",
-	    AdiabaticExpansionInitialVelocity);
-  }
- 
-  return SUCCESS;
-}

File src/enzo/problem_types/CallProblemSpecificRoutines.C

- 
-#ifdef USE_MPI
-#include "communicators.h"
-#endif /* USE_MPI */
- 
-#include "preincludes.h"
-
-#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"
-
-//variables  for RandomForcing
-
-int CallProblemSpecificRoutines(TopGridData * MetaData, HierarchyEntry *ThisGrid,
-				int GridNum, float *norm, float TopGridTimeStep, 
-				int level, int LevelCycleCount[])
-{
-
-  /* Add RandomForcing fields to velocities after the copying of current
-     fields to old. I also update the total energy accordingly here.
-     It makes no sense to force on the very first time step. */
- 
-  if (MetaData->CycleNumber > 0)
-    ThisGrid->GridData->AddRandomForcing(norm, TopGridTimeStep);
-
-  //dcc cut stop Forcing
-
-  if (ProblemType == 24)
-    ThisGrid->GridData->SphericalInfallGetProfile(level, 1);
-  if (ProblemType == 30)
-    ThisGrid->GridData->AnalyzeTrackPeaks(level, 0);
-//  if (ProblemType == 27)
-//    if (ThisGrid->GridData->ReturnProcessorNumber()==MyProcessorNumber){
-//      float AM[3], MeanVelocity[3], DMVelocity[3];
-//      FLOAT Center[] = {0,0,0}, CenterOfMass[3], DMCofM[3];
-//      ThisGrid->GridData->CalculateAngularMomentum
-//	(Center, AM, MeanVelocity, DMVelocity, CenterOfMass, DMCofM);
-//      fprintf(stdout, 
-//	      "level = %"ISYM" %"ISYM" %"ISYM"  "
-//	      "Vel %"FSYM" %"FSYM" %"FSYM"  "
-//	      "DMVel %"FSYM" %"FSYM" %"FSYM"  "
-//	      "CofM %"PSYM" %"PSYM" %"PSYM"  "
-//	      "DMCofM %"FSYM" %"FSYM" %"FSYM"\n",
-//	      level, LevelCycleCount[level], GridNum, MeanVelocity[0],
-//	      MeanVelocity[1], MeanVelocity[2],
-//	      DMVelocity[0], DMVelocity[1], DMVelocity[2],
-//	      -CenterOfMass[0], -CenterOfMass[1], -CenterOfMass[2],
-//	      DMCofM[0], DMCofM[1], DMCofM[2]);
-//    }
-
-  /* Solve analytical free-fall */
-  if (ProblemType == 63) {
-    ThisGrid->GridData->SolveOneZoneFreefall();
-  }
-
-  return SUCCESS;
-}

File src/enzo/problem_types/CheckShearingBoundaryConsistency.C

-/***********************************************************************
-/
-/  CHECK INPUT SHEARING BOUNDARIES ARE CONSISTENT
-/
-/  written by: Fen Zhao
-/  date:       June, 2009
-/  modified1:
-/
-/  PURPOSE:
-/         Check that input parameters for Shearing Boundaries are sensical 
-/         and sets helper variables properlyk 
-/
-************************************************************************/
-
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "TopGridData.h"
-
-int CheckShearingBoundaryConsistency(TopGridData &MetaData){
-
-  //Check only one direction has a shearing boundary condition and both faces shearing
-  int numberOfShearingBoundariesL = 0 ,  numberOfShearingBoundariesR = 0 ;
-  for (int i=0; i<MetaData.TopGridRank; i++){
-    if (MetaData.LeftFaceBoundaryCondition[i]== shearing) { numberOfShearingBoundariesL++; ShearingBoundaryDirection=i;}
-    if (MetaData.RightFaceBoundaryCondition[i]== shearing) { numberOfShearingBoundariesR++; ShearingBoundaryDirection=i;}
-  }
-
-  if (numberOfShearingBoundariesL> 1 || numberOfShearingBoundariesR > 1) 
-    ENZO_FAIL("Too Many Shearing Boundaries");
-  if (numberOfShearingBoundariesL == 0 && numberOfShearingBoundariesR == 0) 
-    return SUCCESS;
- 
-
-  //Check that we have one direction at least where there are periodic BC on both sides, fail otherwise
-  int numberOfPeriodicBoundariesBoth = 0;
-  int tempShearingVelocityDirection;
-  for (int i=0; i<MetaData.TopGridRank; i++)
-      if (MetaData.LeftFaceBoundaryCondition[i] == periodic && MetaData.RightFaceBoundaryCondition[i] == periodic){ 
-      numberOfPeriodicBoundariesBoth++;
-      tempShearingVelocityDirection=i;
-    }
-  
-
-  if (numberOfPeriodicBoundariesBoth == 0) {
-   
-     ENZO_FAIL("Need a Periodic Boundary (on both faces) for a Shearing Boundary");
-  }
-   //if only one direction is periodic, then we autoset ShearingVelocityDirection to the right value
-   if (numberOfPeriodicBoundariesBoth == 1){
-     ShearingVelocityDirection=tempShearingVelocityDirection;
-     return SUCCESS;
-   }
-   
-   //if both directions periodic, then either we use the user specific ShearingVelocityDirection
-   //or use the high index of the two directions
-  
-  if (numberOfPeriodicBoundariesBoth == 2){
-     
-     if (ShearingVelocityDirection==-1) ShearingVelocityDirection=tempShearingVelocityDirection;
-     else if (MetaData.LeftFaceBoundaryCondition[ShearingVelocityDirection]!=periodic 
-	      || MetaData.RightFaceBoundaryCondition[ShearingVelocityDirection]!=periodic)
-       ENZO_FAIL("ShearingVelocityDirection is not Periodic on Both Faces");
-   }
-  
- 
-
-  if (ShearingBoundaryDirection==0){
-    if (ShearingVelocityDirection==1) ShearingOtherDirection=2;
-    else if (ShearingVelocityDirection==2) ShearingOtherDirection=1;}
-  else if (ShearingBoundaryDirection==1){
-    if (ShearingVelocityDirection==0) ShearingOtherDirection=2;
-    else if (ShearingVelocityDirection==2) ShearingOtherDirection=0;}
-  else if (ShearingBoundaryDirection==2){
-    if (ShearingVelocityDirection==0) ShearingOtherDirection=1;
-    else if (ShearingVelocityDirection==1) ShearingOtherDirection=0;} 
-
-
-  return SUCCESS;
-  
-}

File src/enzo/problem_types/ClusterInitialize.C

-/***********************************************************************
-/
-/  INITIALIZE A Cool Core Cluster 
-/
-/  written by: Yuan Li and Greg Bryan
-/  date:       Dec 2011 
-/  modified1:
-/
-/  PURPOSE:
-/
-/  RETURNS: SUCCESS or FAIL
-/
-************************************************************************/
-
-// This routine intializes a new simulation based on the parameter file.
-//
-
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "LevelHierarchy.h"
-#include "TopGridData.h"
-
-void WriteListOfFloats(FILE *fptr, int N, float floats[]);
-void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
-void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
-int RebuildHierarchy(TopGridData *MetaData,
-                     LevelHierarchyEntry *LevelArray[], int level);
-int CommunicationPartitionGrid(HierarchyEntry *Grid, int gridnum);
-
-
-int ClusterInitialize(FILE *fptr, FILE *Outfptr, 
-                           HierarchyEntry &TopGrid, TopGridData &MetaData, ExternalBoundary &Exterior)
-{
-  char *DensName = "Density";
-  char *TEName   = "TotalEnergy";
-  char *GEName   = "GasEnergy";
-  char *Vel1Name = "x-velocity";
-  char *Vel2Name = "y-velocity";
-  char *Vel3Name = "z-velocity";
-  char *ColourName = "colour";
-
-  /* declarations */
-
-  char  line[MAX_LINE_LENGTH];
-  int   dim, ret, level, sphere, i;
-
-  /* set default parameters */
-
-  int ClusterNumberOfSpheres = 1;
-  int ClusterRefineAtStart   = TRUE;
-  int ClusterUseParticles    = FALSE;
-  int ClusterUseColour       = FALSE;
-  float ClusterInitialTemperature = 1000;
-  int   ClusterSphereType[MAX_SPHERES];
-  float ClusterSphereDensity[MAX_SPHERES],
-        ClusterSphereTemperature[MAX_SPHERES],
-        ClusterSphereVelocity[MAX_SPHERES][MAX_DIMENSION],
-        ClusterUniformVelocity[MAX_DIMENSION];
-  FLOAT ClusterSphereRadius[MAX_SPHERES],
-        ClusterSphereCoreRadius[MAX_SPHERES],
-        ClusterSpherePosition[MAX_SPHERES][MAX_DIMENSION];
-
-  for (sphere = 0; sphere < MAX_SPHERES; sphere++) {
-    ClusterSphereRadius[sphere]     = 1.0;
-    ClusterSphereCoreRadius[sphere] = 0.1;
-    ClusterSphereDensity[sphere]    = 1.0;
-    ClusterSphereTemperature[sphere] = 1.0;
-    for (dim = 0; dim < MAX_DIMENSION; dim++) {
-      ClusterSpherePosition[sphere][dim] = 0.5*(DomainLeftEdge[dim] +
-                                                     DomainRightEdge[dim]);
-      ClusterSphereVelocity[sphere][dim] = 0;
-    }
-    ClusterSphereType[sphere]       = 0;
-  }
-  for (dim = 0; dim < MAX_DIMENSION; dim++)
-    ClusterUniformVelocity[dim] = 0;
-
-  /* read input from file */
-
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
-
-    ret = 0;
-
-    /* read parameters */
-
-    ret += sscanf(line, "ClusterNumberOfSpheres = %"ISYM,
-                  &ClusterNumberOfSpheres);
-    ret += sscanf(line, "ClusterRefineAtStart = %"ISYM, 
-                  &ClusterRefineAtStart);
-    ret += sscanf(line, "ClusterUseParticles = %"ISYM, 
-                  &ClusterUseParticles);
-    ret += sscanf(line, "ClusterUseColour = %"ISYM, 
-                  &ClusterUseColour);
-    ret += sscanf(line, "ClusterInitialTemperature = %"FSYM, 
-                  &ClusterInitialTemperature);
-    ret += sscanf(line, "ClusterUniformVelocity = %"FSYM" %"FSYM" %"FSYM, 
-                  ClusterUniformVelocity, ClusterUniformVelocity+1,
-                  ClusterUniformVelocity+2);
-    if (sscanf(line, "ClusterSphereType[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereType[%"ISYM"] = %"ISYM, &sphere,
-                    &ClusterSphereType[sphere]);
-    if (sscanf(line, "ClusterSphereRadius[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereRadius[%"ISYM"] = %"PSYM, &sphere,
-                    &ClusterSphereRadius[sphere]);
-    if (sscanf(line, "ClusterSphereCoreRadius[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereCoreRadius[%"ISYM"] = %"PSYM, &sphere,
-                    &ClusterSphereCoreRadius[sphere]);
-    if (sscanf(line, "ClusterSphereDensity[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereDensity[%"ISYM"] = %"FSYM, &sphere,
-                    &ClusterSphereDensity[sphere]);
-    if (sscanf(line, "ClusterSphereTemperature[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereTemperature[%"ISYM"] = %"FSYM, &sphere,
-                    &ClusterSphereTemperature[sphere]);
-    if (sscanf(line, "ClusterSpherePosition[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSpherePosition[%"ISYM"] = %"PSYM" %"PSYM" %"PSYM, 
-                    &sphere, &ClusterSpherePosition[sphere][0],
-                    &ClusterSpherePosition[sphere][1],
-                    &ClusterSpherePosition[sphere][2]);
-    if (sscanf(line, "ClusterSphereVelocity[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "ClusterSphereVelocity[%"ISYM"] = %"FSYM" %"FSYM" %"FSYM, 
-                    &sphere, &ClusterSphereVelocity[sphere][0],
-                    &ClusterSphereVelocity[sphere][1],
-                    &ClusterSphereVelocity[sphere][2]);
-
-    /* if the line is suspicious, issue a warning */
-
-    if (ret == 0 && strstr(line, "=") && strstr(line, "Cluster") 
-        && line[0] != '#')
-      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
-
-  } // end input from parameter file
-
-  /* set up grid */
-
-  if (TopGrid.GridData->ClusterInitializeGrid(
-             ClusterNumberOfSpheres, ClusterSphereRadius,
-             ClusterSphereCoreRadius, ClusterSphereDensity,
-             ClusterSphereTemperature,
-             ClusterSpherePosition, ClusterSphereVelocity,
-             ClusterSphereType, ClusterUseParticles,
-             ClusterUniformVelocity, ClusterUseColour,
-             ClusterInitialTemperature, 0) == FAIL) {
-    fprintf(stderr, "Error in ClusterInitializeGrid.\n");
-    return FAIL;
-  }
-
-  /* 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]);
-  }
-
-  Exterior.Prepare(TopGrid.GridData);
-
-  float Dummy[MAX_DIMENSION];
-  for (dim = 0; dim < MetaData.TopGridRank; dim++)
-    if (Exterior.InitializeExternalBoundaryFace(dim, 
-                                                MetaData.LeftFaceBoundaryCondition[dim],
-                                                MetaData.RightFaceBoundaryCondition[dim],
-                                                Dummy, Dummy)
-        == FAIL) {
-      fprintf(stderr, "Error in InitializeExternalBoundaryFace.\n");
-      return FAIL;
-    }
-
-  /* Initialize particle boundary conditions. */
-
-  Exterior.InitializeExternalBoundaryParticles(MetaData.ParticleBoundaryType);
-  
-  CommunicationPartitionGrid(&TopGrid, 0);
-
-  /* If requested, refine the grid to the desired level. */
-
-  if (ClusterRefineAtStart) {
-
-    /* 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->ClusterInitializeGrid(
-             ClusterNumberOfSpheres, ClusterSphereRadius,
-             ClusterSphereCoreRadius, ClusterSphereDensity,
-             ClusterSphereTemperature,
-             ClusterSpherePosition, ClusterSphereVelocity,
-             ClusterSphereType, ClusterUseParticles,
-             ClusterUniformVelocity, ClusterUseColour,
-             ClusterInitialTemperature, level+1) == FAIL) {
-          fprintf(stderr, "Error in ClusterInitializeGrid.\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(*Temp->GridHierarchyEntry->ParentGrid->GridData) == FAIL) {
-//                                 *LevelArray[level-1]->GridData) == FAIL) {
-          fprintf(stderr, "Error in grid->ProjectSolutionToParentGrid.\n");
-          return FAIL;
-        }
-        Temp = Temp->NextGridThisLevel;
-      }
-    }
-
-  } // end: if (ClusterRefineAtStart)
-
-  /* set up field names and units */
-
-  int count = 0;
-  DataLabel[count++] = DensName;
-  DataLabel[count++] = TEName;
-  if (DualEnergyFormalism)
-    DataLabel[count++] = GEName;
-  DataLabel[count++] = Vel1Name;
-  DataLabel[count++] = Vel2Name;
-  DataLabel[count++] = Vel3Name;
-  if (ClusterUseColour)
-    DataLabel[count++] = ColourName;
-
-  for (i = 0; i < count; i++)
-    DataUnits[i] = NULL;
-
-  /* Write parameters to parameter output file */
-
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    fprintf(Outfptr, "ClusterNumberOfSpheres    = %"ISYM"\n",
-            ClusterNumberOfSpheres);
-    fprintf(Outfptr, "ClusterRefineAtStart      = %"ISYM"\n",
-            ClusterRefineAtStart);
-    fprintf(Outfptr, "ClusterUseParticles       = %"ISYM"\n",
-            ClusterUseParticles);
-    fprintf(Outfptr, "ClusterUseColour          = %"ISYM"\n",
-            ClusterUseColour);
-    fprintf(Outfptr, "ClusterInitialTemperature = %"GOUTSYM"\n",
-            ClusterInitialTemperature);
-    fprintf(Outfptr, "ClusterUniformVelocity    = %"GOUTSYM" %"GOUTSYM" %"GOUTSYM"\n",
-            ClusterUniformVelocity[0], ClusterUniformVelocity[1],
-            ClusterUniformVelocity[2]);
-    for (sphere = 0; sphere < ClusterNumberOfSpheres; sphere++) {
-      fprintf(Outfptr, "ClusterSphereType[%"ISYM"] = %"ISYM"\n", sphere,
-              ClusterSphereType[sphere]);
-      fprintf(Outfptr, "ClusterSphereRadius[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              ClusterSphereRadius[sphere]);
-      fprintf(Outfptr, "ClusterSphereCoreRadius[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              ClusterSphereCoreRadius[sphere]);
-      fprintf(Outfptr, "ClusterSphereDensity[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              ClusterSphereDensity[sphere]);
-      fprintf(Outfptr, "ClusterSphereTemperature[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              ClusterSphereTemperature[sphere]);
-      fprintf(Outfptr, "ClusterSpherePosition[%"ISYM"] = ", sphere);
-      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
-                        ClusterSpherePosition[sphere]);
-      fprintf(Outfptr, "ClusterSphereVelocity[%"ISYM"] = ", sphere);
-      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
-                        ClusterSphereVelocity[sphere]);
-    }
-  }
-
-  return SUCCESS;
-
-}
-

File src/enzo/problem_types/CollapseTestInitialize.C

-/***********************************************************************
-/
-/  INITIALIZE A COLLAPSE TEST
-/
-/  written by: Greg Bryan
-/  date:       May, 1998
-/  modified1:
-/
-/  PURPOSE:
-/    Set up a number of spherical objects
-/
-/  RETURNS: SUCCESS or FAIL
-/
-************************************************************************/
-
-// This routine intializes a new simulation based on the parameter file.
-//
-
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "LevelHierarchy.h"
-#include "TopGridData.h"
-
-void WriteListOfFloats(FILE *fptr, int N, float floats[]);
-void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
-void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
-int RebuildHierarchy(TopGridData *MetaData,
-		     LevelHierarchyEntry *LevelArray[], int level);
-
-int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, 
-			  HierarchyEntry &TopGrid, TopGridData &MetaData)
-{
-  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 *ColourName = "colour";
-  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 *MetalName = "Metal_Density";
-
-  /* declarations */
-
-  char  line[MAX_LINE_LENGTH];
-  int   dim, ret, level, sphere, i;
-
-  /* set default parameters */
-
-  int CollapseTestNumberOfSpheres = 1;
-  int CollapseTestRefineAtStart   = TRUE;
-  int CollapseTestUseParticles    = FALSE;
-  float CollapseTestParticleMeanDensity = FLOAT_UNDEFINED;
-  int CollapseTestUseColour       = FALSE;
-  int CollapseTestUseMetals       = FALSE;
-  float CollapseTestInitialTemperature = 1000;
-  float CollapseTestInitialDensity     = 1.0;
-  float CollapseTestSphereDensity[MAX_SPHERES],
-    CollapseTestSphereTemperature[MAX_SPHERES],
-    CollapseTestSphereVelocity[MAX_SPHERES][MAX_DIMENSION],
-    CollapseTestUniformVelocity[MAX_DIMENSION],
-    CollapseTestFracKeplerianRot[MAX_SPHERES],
-    CollapseTestSphereTurbulence[MAX_SPHERES],
-    CollapseTestSphereDispersion[MAX_SPHERES],
-    CollapseTestSphereCutOff[MAX_SPHERES],
-    CollapseTestSphereAng1[MAX_SPHERES],
-    CollapseTestSphereAng2[MAX_SPHERES],
-    CollapseTestSphereRotationPeriod[MAX_SPHERES],
-    CollapseTestSphereMetallicity[MAX_SPHERES];
-  int CollapseTestSphereNumShells[MAX_SPHERES],
-    CollapseTestSphereInitialLevel[MAX_SPHERES],
-    CollapseTestSphereType[MAX_SPHERES];
-  FLOAT CollapseTestSphereRadius[MAX_SPHERES],
-    CollapseTestSphereCoreRadius[MAX_SPHERES],
-    CollapseTestSpherePosition[MAX_SPHERES][MAX_DIMENSION];
-
-  for (sphere = 0; sphere < MAX_SPHERES; sphere++) {
-    CollapseTestSphereRadius[sphere]     = 1.0;
-    CollapseTestSphereCoreRadius[sphere] = 0.1;
-    CollapseTestSphereDensity[sphere]    = 1.0;
-    CollapseTestSphereTemperature[sphere] = 1.0;
-    CollapseTestFracKeplerianRot[sphere] = 0.0;
-    CollapseTestSphereTurbulence[sphere] = 0.0;
-    CollapseTestSphereDispersion[sphere] = 0.0;
-    CollapseTestSphereCutOff[sphere] = 6.5;
-    CollapseTestSphereAng1[sphere] = 0;
-    CollapseTestSphereAng2[sphere] = 0;
-    CollapseTestSphereRotationPeriod[sphere] = 0;
-    CollapseTestSphereNumShells[sphere] = 1;
-    CollapseTestSphereMetallicity[sphere] = 0;
-    CollapseTestSphereInitialLevel[sphere] = 0;
-
-    for (dim = 0; dim < MAX_DIMENSION; dim++) {
-      CollapseTestSpherePosition[sphere][dim] = 0.5*(DomainLeftEdge[dim] +
-						     DomainRightEdge[dim]);
-      CollapseTestSphereVelocity[sphere][dim] = 0;
-    }
-    CollapseTestSphereType[sphere]       = 0;
-  }
-  for (dim = 0; dim < MAX_DIMENSION; dim++)
-    CollapseTestUniformVelocity[dim] = 0;
-
-  /* read input from file */
-
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
-
-    ret = 0;
-
-    /* read parameters */
-
-    ret += sscanf(line, "CollapseTestNumberOfSpheres = %"ISYM,
-		  &CollapseTestNumberOfSpheres);
-    ret += sscanf(line, "CollapseTestRefineAtStart = %"ISYM, 
-		  &CollapseTestRefineAtStart);
-    ret += sscanf(line, "CollapseTestUseParticles = %"ISYM, 
-		  &CollapseTestUseParticles);
-    ret += sscanf(line, "CollapseTestParticleMeanDensity = %"FSYM,
-		  &CollapseTestParticleMeanDensity);
-    ret += sscanf(line, "CollapseTestUseColour = %"ISYM, 
-		  &CollapseTestUseColour);
-    ret += sscanf(line, "CollapseTestUseMetals = %"ISYM, 
-		  &CollapseTestUseMetals);
-    ret += sscanf(line, "CollapseTestInitialTemperature = %"FSYM, 
-		  &CollapseTestInitialTemperature);
-    ret += sscanf(line, "CollapseTestInitialDensity = %"FSYM,
-		  &CollapseTestInitialDensity);
-    ret += sscanf(line, "CollapseTestUniformVelocity = %"FSYM" %"FSYM" %"FSYM, 
-		  CollapseTestUniformVelocity, CollapseTestUniformVelocity+1,
-		  CollapseTestUniformVelocity+2);
-    if (sscanf(line, "CollapseTestSphereType[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereType[%"ISYM"] = %"ISYM, &sphere,
-		    &CollapseTestSphereType[sphere]);
-    if (sscanf(line, "CollapseTestSphereRadius[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereRadius[%"ISYM"] = %"PSYM, &sphere,
-		    &CollapseTestSphereRadius[sphere]);
-    if (sscanf(line, "CollapseTestSphereCoreRadius[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereCoreRadius[%"ISYM"] = %"PSYM, &sphere,
-		    &CollapseTestSphereCoreRadius[sphere]);
-    if (sscanf(line, "CollapseTestSphereDensity[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereDensity[%"ISYM"] = %"FSYM, &sphere,
-		    &CollapseTestSphereDensity[sphere]);
-    if (sscanf(line, "CollapseTestSphereTemperature[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereTemperature[%"ISYM"] = %"FSYM, &sphere,
-		    &CollapseTestSphereTemperature[sphere]);
-    if (sscanf(line, "CollapseTestSphereMetallicity[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereMetallicity[%"ISYM"] = %"FSYM, &sphere,
-		    &CollapseTestSphereMetallicity[sphere]);
-    if (sscanf(line, "CollapseTestSpherePosition[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSpherePosition[%"ISYM"] = %"PSYM" %"PSYM" %"PSYM, 
-		    &sphere, &CollapseTestSpherePosition[sphere][0],
-		    &CollapseTestSpherePosition[sphere][1],
-		    &CollapseTestSpherePosition[sphere][2]);
-    if (sscanf(line, "CollapseTestSphereVelocity[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereVelocity[%"ISYM"] = %"FSYM" %"FSYM" %"FSYM, 
-		    &sphere, &CollapseTestSphereVelocity[sphere][0],
-		    &CollapseTestSphereVelocity[sphere][1],
-		    &CollapseTestSphereVelocity[sphere][2]);
-    if (sscanf(line, "CollapseTestFracKeplerianRot[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestFracKeplerianRot[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestFracKeplerianRot[sphere]);
-    if (sscanf(line, "CollapseTestSphereTurbulence[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereTurbulence[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestSphereTurbulence[sphere]);
-    if (sscanf(line, "CollapseTestSphereDispersion[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereDispersion[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestSphereDispersion[sphere]);
-    if (sscanf(line, "CollapseTestSphereCutOff[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereCutOff[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestSphereCutOff[sphere]);
-    if (sscanf(line, "CollapseTestSphereAng1[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereAng1[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestSphereAng1[sphere]);
-    if (sscanf(line, "CollapseTestSphereAng2[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereAng2[%"ISYM"] = %"FSYM, &sphere,
-                    &CollapseTestSphereAng2[sphere]);
-    if (sscanf(line, "CollapseTestSphereRotationPeriod[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereRotationPeriod[%"ISYM"] = %"FSYM, &sphere,
-		    &CollapseTestSphereRotationPeriod[sphere]);
-    if (sscanf(line, "CollapseTestSphereNumShells[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereNumShells[%"ISYM"] = %"ISYM, &sphere,
-                    &CollapseTestSphereNumShells[sphere]);
-    if (sscanf(line, "CollapseTestSphereInitialLevel[%"ISYM"]", &sphere) > 0)
-      ret += sscanf(line, "CollapseTestSphereInitialLevel[%"ISYM"] = %"ISYM, &sphere,
-                    &CollapseTestSphereInitialLevel[sphere]);
-
-    /* if the line is suspicious, issue a warning */
-
-    if (ret == 0 && strstr(line, "=") && strstr(line, "CollapseTest") 
-	&& line[0] != '#')
-      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
-
-  } // end input from parameter file
-
-  /* set up grid */
-
-  if (TopGrid.GridData->CollapseTestInitializeGrid(
-	     CollapseTestNumberOfSpheres, CollapseTestSphereRadius,
-	     CollapseTestSphereCoreRadius, CollapseTestSphereDensity,
-	     CollapseTestSphereTemperature, CollapseTestSphereMetallicity,
-	     CollapseTestSpherePosition, CollapseTestSphereVelocity,
-             CollapseTestFracKeplerianRot, CollapseTestSphereTurbulence,
-	     CollapseTestSphereDispersion,
-             CollapseTestSphereCutOff, CollapseTestSphereAng1,
-             CollapseTestSphereAng2, CollapseTestSphereRotationPeriod,
-	     CollapseTestSphereNumShells,
-	     CollapseTestSphereType, CollapseTestUseParticles,
-	     CollapseTestParticleMeanDensity,
-             CollapseTestUniformVelocity, CollapseTestUseColour,
-	     CollapseTestUseMetals,
-             CollapseTestInitialTemperature, CollapseTestInitialDensity,
-	     0) == FAIL) {
-    ENZO_FAIL("Error in CollapseTestInitializeGrid.");
-  }
-
-  /* Convert minimum initial overdensity for refinement to mass
-     (unless MinimumMass itself was actually set). */
-
-  for (i = 0; i < MAX_FLAGGING_METHODS; i++)
-    if (MinimumMassForRefinement[i] == FLOAT_UNDEFINED) {
-      MinimumMassForRefinement[i] = MinimumOverDensityForRefinement[i];
-      for (dim = 0; dim < MetaData.TopGridRank; dim++)
-	MinimumMassForRefinement[i] *=
-	  (DomainRightEdge[dim]-DomainLeftEdge[dim])/
-	  float(MetaData.TopGridDims[dim]);
-    }
-
-  /* If requested and there are no manual settings of the refinement
-     of spheres, refine the grid to the desired level. */
-
-  int MaxInitialLevel = 0;
-  for (sphere = 0; sphere < CollapseTestNumberOfSpheres; sphere++)
-    MaxInitialLevel = max(MaxInitialLevel, CollapseTestSphereInitialLevel[sphere]);
-
-  if (CollapseTestRefineAtStart) {
-
-    /* If the user specified an initial refinement level for a sphere,
-       then manually create the hierarchy first. */
-
-    if (MaxInitialLevel > 0) {
-
-      int lev, max_level;
-      float dx;
-      HierarchyEntry **Subgrid;
-      int NumberOfSubgridDims[MAX_DIMENSION];
-      FLOAT ThisLeftEdge[MAX_DIMENSION], ThisRightEdge[MAX_DIMENSION];
-
-      for (sphere = 0; sphere < CollapseTestNumberOfSpheres; sphere++) {
-	
-	max_level = CollapseTestSphereInitialLevel[sphere];
-	if (max_level > 0) {
-
-	  Subgrid = new HierarchyEntry*[max_level];
-	  for (lev = 0; lev < max_level; lev++)
-	    Subgrid[lev] = new HierarchyEntry;
-
-	  for (lev = 0; lev < max_level; lev++) {
-	    
-	    for (dim = 0; dim < MetaData.TopGridRank; dim++) {
-	      dx = 1.0 / float(MetaData.TopGridDims[dim]) / POW(RefineBy, lev);
-	      ThisLeftEdge[dim] = CollapseTestSpherePosition[sphere][dim] -
-		0.5 * CollapseTestSphereRadius[sphere] - 2*dx;  // plus some buffer
-	      ThisLeftEdge[dim] = nint(ThisLeftEdge[dim] / dx) * dx;
-	      ThisRightEdge[dim] = CollapseTestSpherePosition[sphere][dim] +
-		0.5 * CollapseTestSphereRadius[sphere] + 2*dx;
-	      ThisRightEdge[dim] = nint(ThisRightEdge[dim] / dx) * dx;
-	      NumberOfSubgridDims[dim] = 
-		nint((ThisRightEdge[dim] - ThisLeftEdge[dim]) / 
-		     (DomainRightEdge[dim] - DomainLeftEdge[dim]) / dx);		
-	    } // ENDFOR dims
-
-	    if (debug)
-	      printf("CollapseTest:: Level[%"ISYM"]: NumberOfSubgridZones[0] = %"ISYM"\n",
-		     lev+1, NumberOfSubgridDims[0]);
-	    
-	    if (NumberOfSubgridDims[0] > 0) {
-
-	      // Insert into AMR hierarchy
-	      if (lev == 0) {
-		Subgrid[lev]->NextGridThisLevel = TopGrid.NextGridNextLevel;
-		TopGrid.NextGridNextLevel = Subgrid[lev];
-		Subgrid[lev]->ParentGrid = &TopGrid;
-	      } else {
-		Subgrid[lev]->NextGridThisLevel = NULL;
-		Subgrid[lev]->ParentGrid = Subgrid[lev-1];
-	      }
-	      if (lev == max_level-1)
-		Subgrid[lev]->NextGridNextLevel = NULL;
-	      else
-		Subgrid[lev]->NextGridNextLevel = Subgrid[lev+1];
-
-	      // Create grid
-	      for (dim = 0; dim < MetaData.TopGridRank; dim++)
-		NumberOfSubgridDims[dim] += 2*DEFAULT_GHOST_ZONES;
-	      Subgrid[lev]->GridData = new grid;
-	      Subgrid[lev]->GridData->InheritProperties(TopGrid.GridData);
-	      Subgrid[lev]->GridData->PrepareGrid(MetaData.TopGridRank, 
-						  NumberOfSubgridDims,
-						  ThisLeftEdge,
-						  ThisRightEdge, 0);
-
-
-	      if (Subgrid[lev]->GridData->CollapseTestInitializeGrid(
-	          CollapseTestNumberOfSpheres, CollapseTestSphereRadius,
-		  CollapseTestSphereCoreRadius, CollapseTestSphereDensity,
-		  CollapseTestSphereTemperature, CollapseTestSphereMetallicity,
-		  CollapseTestSpherePosition, CollapseTestSphereVelocity,
-		  CollapseTestFracKeplerianRot, CollapseTestSphereTurbulence,
-		  CollapseTestSphereDispersion,
-		  CollapseTestSphereCutOff, CollapseTestSphereAng1,
-		  CollapseTestSphereAng2, CollapseTestSphereRotationPeriod,
-		  CollapseTestSphereNumShells,
-		  CollapseTestSphereType, CollapseTestUseParticles,
-		  CollapseTestParticleMeanDensity,
-		  CollapseTestUniformVelocity, CollapseTestUseColour,
-		  CollapseTestUseMetals,
-		  CollapseTestInitialTemperature, CollapseTestInitialDensity,
-		  lev-1) == FAIL) {
-		ENZO_FAIL("Error in CollapseTestInitializeGrid.");
-	      }
-	      
-	    } // ENDIF zones exist
-	  } // ENDFOR levels
-	} // ENDIF max_level > 0
-      } // ENDFOR spheres
-    } // ENDIF MaxInitialLevel > 0
-
-    /* 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. */
-
-    if (MaxInitialLevel == 0) {
-      for (level = 0; level < MaximumRefinementLevel; level++) {
-	if (RebuildHierarchy(&MetaData, LevelArray, level) == FAIL) {
-	  ENZO_FAIL("Error in RebuildHierarchy.");
-	}
-	if (LevelArray[level+1] == NULL)
-	  break;
-	LevelHierarchyEntry *Temp = LevelArray[level+1];
-	while (Temp != NULL) {
-	  if (Temp->GridData->CollapseTestInitializeGrid(
-		 CollapseTestNumberOfSpheres, CollapseTestSphereRadius,
-		 CollapseTestSphereCoreRadius, CollapseTestSphereDensity,
-		 CollapseTestSphereTemperature, CollapseTestSphereMetallicity,
-		 CollapseTestSpherePosition, CollapseTestSphereVelocity,
-		 CollapseTestFracKeplerianRot, CollapseTestSphereTurbulence,
-		 CollapseTestSphereDispersion,
-		 CollapseTestSphereCutOff, CollapseTestSphereAng1,
-		 CollapseTestSphereAng2, CollapseTestSphereRotationPeriod,
-		 CollapseTestSphereNumShells,
-		 CollapseTestSphereType, CollapseTestUseParticles,
-		 CollapseTestParticleMeanDensity,
-		 CollapseTestUniformVelocity, CollapseTestUseColour,
-		 CollapseTestUseMetals,
-		 CollapseTestInitialTemperature, CollapseTestInitialDensity,
-		 level+1) == FAIL) {
-	    ENZO_FAIL("Error in CollapseTestInitializeGrid.");
-	  }
-	  Temp = Temp->NextGridThisLevel;
-	}
-      } // end: loop over levels
-    } // ENDELSE manually set refinement 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) {
-	  ENZO_FAIL("Error in grid->ProjectSolutionToParentGrid.");
-	}
-	Temp = Temp->NextGridThisLevel;
-      }
-    }
-
-  } // end: if (CollapseTestRefineAtStart)
-
-  /* set up field names and units */
-
-  int count = 0;
-  DataLabel[count++] = (char*) DensName;
-  DataLabel[count++] = (char*) TEName;
-  if (DualEnergyFormalism)
-    DataLabel[count++] = (char*) GEName;
-  DataLabel[count++] = (char*) Vel1Name;
-  DataLabel[count++] = (char*) Vel2Name;
-  DataLabel[count++] = (char*) Vel3Name;
-  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 (CollapseTestUseColour)
-    DataLabel[count++] = (char*) ColourName;
-  if (CollapseTestUseMetals)
-    DataLabel[count++] = (char*) MetalName;
-
-  for (i = 0; i < count; i++)
-    DataUnits[i] = NULL;
-
-  /* Write parameters to parameter output file */
-
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    fprintf(Outfptr, "CollapseTestNumberOfSpheres    = %"ISYM"\n",
-	    CollapseTestNumberOfSpheres);
-    fprintf(Outfptr, "CollapseTestRefineAtStart      = %"ISYM"\n",
-	    CollapseTestRefineAtStart);
-    fprintf(Outfptr, "CollapseTestUseParticles       = %"ISYM"\n",
-	    CollapseTestUseParticles);
-    fprintf(Outfptr, "CollapseTestUseColour          = %"ISYM"\n",
-	    CollapseTestUseColour);
-    fprintf(Outfptr, "CollapseTestUseMetals          = %"ISYM"\n",
-	    CollapseTestUseMetals);
-    fprintf(Outfptr, "CollapseTestInitialTemperature = %"FSYM"\n",
-	    CollapseTestInitialTemperature);
-    fprintf(Outfptr, "CollapseTestInitialDensity     = %"FSYM"\n",
-	    CollapseTestInitialDensity);
-    fprintf(Outfptr, "CollapseTestUniformVelocity    = %"FSYM" %"FSYM" %"FSYM"\n",
-	    CollapseTestUniformVelocity[0], CollapseTestUniformVelocity[1],
-	    CollapseTestUniformVelocity[2]);
-    for (sphere = 0; sphere < CollapseTestNumberOfSpheres; sphere++) {
-      fprintf(Outfptr, "CollapseTestSphereType[%"ISYM"] = %"ISYM"\n", sphere,
-	      CollapseTestSphereType[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereRadius[%"ISYM"] = %"GOUTSYM"\n", sphere,
-	      CollapseTestSphereRadius[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereCoreRadius[%"ISYM"] = %"GOUTSYM"\n", sphere,
-	      CollapseTestSphereCoreRadius[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereDensity[%"ISYM"] = %"FSYM"\n", sphere,
-	      CollapseTestSphereDensity[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereTemperature[%"ISYM"] = %"FSYM"\n", sphere,
-	      CollapseTestSphereTemperature[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereMetallicity[%"ISYM"] = %"FSYM"\n", sphere,
-	      CollapseTestSphereMetallicity[sphere]);
-      fprintf(Outfptr, "CollapseTestSpherePosition[%"ISYM"] = ", sphere);
-      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
-			CollapseTestSpherePosition[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereVelocity[%"ISYM"] = ", sphere);
-      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
-			CollapseTestSphereVelocity[sphere]);
-      fprintf(Outfptr, "CollapseTestFracKeplerianRot[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              CollapseTestFracKeplerianRot[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereTurbulence[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              CollapseTestSphereTurbulence[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereCutOff[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              CollapseTestSphereCutOff[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereAng1[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              CollapseTestSphereAng1[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereAng2[%"ISYM"] = %"GOUTSYM"\n", sphere,
-              CollapseTestSphereAng2[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereRotationPeriod[%"ISYM"] = %"GOUTSYM"\n", sphere,
-	      CollapseTestSphereRotationPeriod[sphere]);
-      fprintf(Outfptr, "CollapseTestSphereNumShells[%"ISYM"] = %"ISYM"\n", sphere,
-              CollapseTestSphereNumShells[sphere]);
-    }
-  }
-
-  return SUCCESS;
-
-}

File src/enzo/problem_types/ComputeRandomForcingNormalization.C

-/***********************************************************************
-/
-/  COMPUTE RANDOM FORCING NORMALIZATION
-/
-/  written by: Alexei Kritsuk
-/  date:       January, 2004
-/  modified1:
-/
-/  PURPOSE:
-/
-************************************************************************/
-#ifdef USE_MPI
-#include "communicators.h"
-#endif /* USE_MPI */
-
-#include "preincludes.h"
-#include "ErrorExceptions.h"
-#include "performance.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"
- 
-/* ======================================================================= */
-/* Function prototypes. */
- 
-int   GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level,
-			HierarchyEntry **Grids[]);
- 
-/* This routine calculates the normalization for Random Forcing on
-   level 0 (since this involves communication). */
- 
-int ComputeRandomForcingNormalization(LevelHierarchyEntry *LevelArray[],
-				      int level, TopGridData *MetaData,
-				      float * norm, float * pTopGridTimeStep)
-{
-  /* Return if this does not concern us */
-  if (!RandomForcing) return SUCCESS;
- 
-  LCAPERF_START("ComputeRandomForcingNormalization");
-
-  /* If level is above 0 then complain: forcing will only work on level 0
-     grid(s). */
-
-  if ((MetaData->CycleNumber <= 0) || (level != 0)) {
-      LCAPERF_STOP("ComputeRandomForcingNormalization");
-      return SUCCESS;
-  }
- 
-  /* Create an array (Grids) of all the grids on level 0. */
- 
-  typedef HierarchyEntry* HierarchyEntryPointer;
-  HierarchyEntry **Grids;
-  int NumberOfGrids = GenerateGridArray(LevelArray, level, &Grids);
- 
-  /* Loop over level 0 grids and compute sums over each of the grids. */
- 
-  int grid, GlobNum = 9;
-  float * GlobVal = new float[GlobNum];
-  for (int num = 0; num < GlobNum; num++)
-    GlobVal[num] = 0.0;
-  for (grid = 0; grid < NumberOfGrids; grid++)
-    if (Grids[grid]->GridData->PrepareRandomForcingNormalization(GlobVal,
-								 GlobNum)
-	== FAIL) {
-      ENZO_FAIL("Error in grid->PrepareRandomForcingNormalization.");
-    }
- 
-  /* Communicate grid-specific sums and compute global sums;
-     also communicate min/max Density values. */
- 
-  CommunicationAllSumValues(GlobVal, GlobNum-2);
-  float minDens = CommunicationMinValue(*(GlobVal+GlobNum-2));
-  float maxDens = CommunicationMaxValue(*(GlobVal+GlobNum-1));
- 
-  /* Compute normalization. */
- 
-  int numberOfGridZones = 1;
-  for (int dim = 0; dim <  MetaData->TopGridRank; dim++)
-    numberOfGridZones *= MetaData->TopGridDims[dim];
- 
-  float dt = LevelArray[level]->GridData->ReturnTimeStep();
-  *pTopGridTimeStep = dt;
-  if (RandomForcingEdot == 0.0)
-    *norm = 0.0;
-  else
-    *norm = ( sqrt(GlobVal[0]*GlobVal[0] + GlobVal[1]*dt*RandomForcingEdot*2.0*
-		 numberOfGridZones) - GlobVal[0] )/GlobVal[1];
- 
-  if (debug) printf("RandomForcingNormalization %.10"GSYM"\n", *norm);
-  if (debug) printf("RandomForcingGlobals: E_k %"GSYM" M_m %"GSYM" M_v %"GSYM" \n",
-		    0.5*GlobVal[4]/numberOfGridZones,
-		    sqrt(GlobVal[2]/numberOfGridZones),
-		    sqrt(GlobVal[3]/numberOfGridZones));
- 
- 
-  /* Output results for level 0 grids. */
- 
-  if (MetaData->CycleSkipGlobalDataDump != 0)
-    if (MyProcessorNumber == ROOT_PROCESSOR &&
-	MetaData->CycleNumber % MetaData->CycleSkipGlobalDataDump == 0.0) {
-      FILE * Fptr;
-      if ((Fptr = fopen("randomForcing.out", "a")) == NULL)
-	ERROR_MESSAGE;
-      fprintf( Fptr, "%"ISYM" %9.6"FSYM" ", MetaData->CycleNumber, MetaData->Time);
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[1]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[2]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[3]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[4]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[5]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[6]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[7]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[8]); //gv1
-      fprintf( Fptr, " %9.3"GSYM, GlobVal[9]); //gv1
-      fprintf( Fptr, " %9.6"FSYM, 0.50*GlobVal[4]/numberOfGridZones);   // kinetic energy
-      fprintf( Fptr, " %9.6"FSYM, sqrt(GlobVal[2]/numberOfGridZones));  // mass weighted rms Mach
-      fprintf( Fptr, " %9.6"FSYM, sqrt(GlobVal[3]/numberOfGridZones));  // volume weighed rms Mach
-      fprintf( Fptr, " %9.6"FSYM, sqrt(GlobVal[5]/numberOfGridZones));  // rms Velocity
-      fprintf( Fptr, " %9.6"FSYM, sqrt(GlobVal[6]/numberOfGridZones));  // Density variance
-      fprintf( Fptr, " %9.6"FSYM" %10.5"FSYM"\n", minDens, maxDens );                  // min/max Density
-      fclose(Fptr);
-    }
- 
-  /* clean up. */
- 
-  delete [] GlobVal;
- 
-  LCAPERF_STOP("ComputeRandomForcingNormalization");
-  return SUCCESS;
-}

File src/enzo/problem_types/ConductionBubbleInitialize.C

-////////////////////////////////////////////////////////////////////////////////
-//
-//  Conduction Test Problem
-//
-//  written by: David A. Ventimiglia, Brian O'Shea
-//  date:       June 2009
-//  modified:  February 10, 2011 by BWO
-//
-//  PURPOSE: This initializes a bubble with user-controlled entropy in an 
-//     ambient medium that is initially in hydrostatic equilibrium
-//     (with density, temperature profile controlled by the user).  The
-//     bubble then does its thing based on its entropy relative to the 
-//     ambient medium.
-//
-//  RETURNS: SUCCESS or FAIL
-//
-////////////////////////////////////////////////////////////////////////////////
- 
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, double *MassUnits, FLOAT Time);
-
-// Problem Initializer
-int ConductionBubbleInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData){
-
-  if(debug){
-    printf("Entering ConductionBubbleInitialize\n");
-    fflush(stdout);
-  }
-
-  char line[MAX_LINE_LENGTH];
-  float LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
-  int i, j, dim, ret;
-
-  float ConductionBubbleDensity = 1.0;
-  float ConductionBubbleTotalEnergy = 1.0;
-  float ConductionBubbleGasEnergy = 1.0;
-  float ConductionBubbleVelocity[3] = {0.0,0.0,0.0};
-  float ConductionBubbleInitialUniformBField[3] = {0.0,0.0,0.0};  // in Gauss
-
-  FLOAT ConductionBubbleRadiusOfBubble = 0.1;  // units of box size
-  int   ConductionBubblePulseType = 1;  // pulse type
-  float ConductionBubbleDeltaEntropy = 0.1;    // multiple of entropy
-  float ConductionBubbleMidpointEntropy = 20.0;  // entropy in kev*cm^2
-  float ConductionBubbleEntropyGradient = 1.0;   // kev*cm^2 / kpc
-  float ConductionBubbleMidpointTemperature = 5.0e+7;  // Kelvin
-  FLOAT ConductionBubbleCenter[MAX_DIMENSION] = {0.5,0.5,0.5};
-
-  // Read parameters
-
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
-    ret = 0;
-    ret += sscanf(line, "ConductionBubbleRadiusOfBubble = %"PSYM, &ConductionBubbleRadiusOfBubble);
-    ret += sscanf(line, "ConductionBubblePulseType = %"ISYM, &ConductionBubblePulseType);
-    ret += sscanf(line, "ConductionBubbleDeltaEntropy = %"FSYM, &ConductionBubbleDeltaEntropy);
-    ret += sscanf(line, "ConductionBubbleMidpointEntropy = %"FSYM, &ConductionBubbleMidpointEntropy);
-    ret += sscanf(line, "ConductionBubbleEntropyGradient = %"FSYM, &ConductionBubbleEntropyGradient);
-    ret += sscanf(line, "ConductionBubbleMidpointTemperature = %"FSYM, &ConductionBubbleMidpointTemperature);
-    ret += sscanf(line, "ConductionBubbleCenter = %"PSYM" %"PSYM" %"PSYM, &ConductionBubbleCenter[0],
-		  &ConductionBubbleCenter[1],&ConductionBubbleCenter[2]);
-    ret += sscanf(line, "TestProblemUseMetallicityField  = %"ISYM, &TestProblemData.UseMetallicityField);
-    ret += sscanf(line, "TestProblemInitialMetallicityFraction  = %"FSYM, &TestProblemData.MetallicityField_Fraction);
-    ret += sscanf(line, "ConductionBubbleBField = %"FSYM" %"FSYM" %"FSYM,&ConductionBubbleInitialUniformBField[0],
-		  &ConductionBubbleInitialUniformBField[1], &ConductionBubbleInitialUniformBField[2]);
-
-    if (ret == 0 && 
-	strstr(line, "=") && strstr(line, "ConductionBubble") &&
-	line[0] != '#' && MyProcessorNumber == ROOT_PROCESSOR) {
-      fprintf(stderr, "*** warning: the following parameter line was not interpreted:\n%s\n", line);
-    }
-  }
-
-  /* error checking */
-  if (Mu != 0.6) {
-    if (MyProcessorNumber == ROOT_PROCESSOR)
-      fprintf(stderr, "warning: mu = 0.6 assumed in initialization; setting Mu = 0.6 for consistency.\n");
-    Mu = 0.6;
-  }
-
-  ConductionBubbleGasEnergy = ConductionBubbleTotalEnergy;
-
-  float DensityUnits=1.0, LengthUnits=1.0, TemperatureUnits=1.0, TimeUnits=1.0,
-    VelocityUnits=1.0;
-  double MassUnits=1.0;
-
-  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	       &TimeUnits, &VelocityUnits, &MassUnits, 0.0) == FAIL) {
-    fprintf(stderr, "Error in GetUnits.\n");
-    return FAIL;
-  }
-
-  /* if we're using MHD, we assume the user inputs values of the B-field.
-     Convert to internal units, and add the magnetic field energy to the 
-     total energy (that's not done automatically, since we don't know for sure
-     what the user is going to do). */
-  if (HydroMethod == MHD_RK){
-    float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
-    ConductionBubbleInitialUniformBField[0] /= MagneticUnits;
-    ConductionBubbleInitialUniformBField[1] /= MagneticUnits;
-    ConductionBubbleInitialUniformBField[2] /= MagneticUnits;
-
-    ConductionBubbleTotalEnergy += 0.5*(ConductionBubbleInitialUniformBField[0]*ConductionBubbleInitialUniformBField[0] + 
-				      ConductionBubbleInitialUniformBField[1]*ConductionBubbleInitialUniformBField[1] + 
-				      ConductionBubbleInitialUniformBField[2]*ConductionBubbleInitialUniformBField[2])/ConductionBubbleDensity;
-  }
-
-  // Create a uniform grid
-  if (TopGrid.GridData->InitializeUniformGrid(ConductionBubbleDensity,
-					      ConductionBubbleTotalEnergy,
-					      ConductionBubbleGasEnergy,
-					      ConductionBubbleVelocity,
-					      ConductionBubbleInitialUniformBField) == FAIL) {
-    ENZO_FAIL("Error in InitializeUniformGrid.");
-  }
-  
-  // Then perturb it
-  if (TopGrid.GridData->ConductionBubbleInitialize(ConductionBubbleRadiusOfBubble, 
-						   ConductionBubblePulseType,
-						   ConductionBubbleDeltaEntropy, 
-						   ConductionBubbleMidpointEntropy, 
-						   ConductionBubbleEntropyGradient, 
-						   ConductionBubbleMidpointTemperature,
-						   ConductionBubbleCenter) == FAIL) {
-    ENZO_FAIL("Error in ConductionBubbleInitialize.");
-  }
-
-  // set up field names and units
-  i = 0;
-  DataLabel[i++] = "Density";
-  DataLabel[i++] = "TotalEnergy";
-  if (DualEnergyFormalism) {DataLabel[i++] = "GasEnergy";}
-  if (MetaData.TopGridRank > 0) {DataLabel[i++] = "x-velocity";}
-  if (MetaData.TopGridRank > 1 || HydroMethod > 2) {DataLabel[i++] = "y-velocity";}
-  if (MetaData.TopGridRank > 2 || HydroMethod > 2) {DataLabel[i++] = "z-velocity";}
-  if (HydroMethod == MHD_RK) {
-    DataLabel[i++] = "Bx";
-    DataLabel[i++] = "By";
-    DataLabel[i++] = "Bz";
-    DataLabel[i++] = "Phi";
-    if(UseDivergenceCleaning){
-      DataLabel[i++] = "Phip";
-    }
-  }
-
-  if (TestProblemData.UseMetallicityField)
-    DataLabel[i++] = "Metal_Density";
-
-  for (j=0; j < i; j++) 
-    DataUnits[j] = NULL;
-
-  /* Write parameters to parameter output file */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    fprintf(Outfptr, "ConductionBubbleRadiusOfBubble = %"PSYM"\n", ConductionBubbleRadiusOfBubble);
-    fprintf(Outfptr, "ConductionBubblePulseType = %"ISYM"\n", ConductionBubblePulseType);
-    fprintf(Outfptr, "ConductionBubbleDeltaEntropy = %"FSYM"\n", ConductionBubbleDeltaEntropy);
-    fprintf(Outfptr, "ConductionBubbleMidpointEntropy = %"FSYM"\n", ConductionBubbleMidpointEntropy);
-    fprintf(Outfptr, "ConductionBubbleEntropyGradient = %"FSYM"\n", ConductionBubbleEntropyGradient);
-    fprintf(Outfptr, "ConductionBubbleMidpointTemperature = %"FSYM"\n", ConductionBubbleMidpointTemperature);
-    fprintf(Outfptr, "ConductionBubbleCenter = %"PSYM" %"PSYM" %"PSYM"\n", ConductionBubbleCenter,
-		  ConductionBubbleCenter+1,ConductionBubbleCenter+2);
-    fprintf(Outfptr, "TestProblemUseMetallicityField  = %"ISYM"\n", TestProblemData.UseMetallicityField);
-    fprintf(Outfptr, "TestProblemInitialMetallicityFraction  = %"FSYM"\n", TestProblemData.MetallicityField_Fraction);
-    fprintf(Outfptr, "ConductionBubbleBField = %"FSYM" %"FSYM" %"FSYM"\n",ConductionBubbleInitialUniformBField[0],
-	    ConductionBubbleInitialUniformBField[1],ConductionBubbleInitialUniformBField[2]);
-  }
-
-  if(debug){
-    printf("Exiting ConductionBubbleInitialize\n");
-    fflush(stdout);
-  }
-
-  return SUCCESS;
-
-}

File src/enzo/problem_types/ConductionCloudInitialize.C

-////////////////////////////////////////////////////////////////////////////////
-//
-//  Conduction Test Problem
-//
-//  written by: David A. Ventimiglia, Brian O'Shea
-//  date:       June 2009
-//  modified:  
-//
-//  PURPOSE: 
-//
-//  RETURNS: SUCCESS or FAIL
-//
-////////////////////////////////////////////////////////////////////////////////
- 
-#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 "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-
-int GetUnits(float *DensityUnits, float *LengthUnits,
-	     float *TemperatureUnits, float *TimeUnits,
-	     float *VelocityUnits, double *MassUnits, FLOAT Time);
-
-// Problem Initializer
-int ConductionCloudInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, TopGridData &MetaData){
-
-  if(debug){
-    printf("Entering ConductionCloudInitialize\n");
-    fflush(stdout);
-  }
-
-  char line[MAX_LINE_LENGTH];
-  float LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
-  int i, j, dim, ret;
-
-  float ConductionCloudDensity = 1.0;
-  float ConductionCloudTemperature = 1.0;
-  float ConductionCloudTotalEnergy = 1.0;
-  float ConductionCloudVelocity[3] = {0.0,0.0,0.0};
-  float ConductionCloudInitialUniformBField[3] = {0.0,0.0,0.0};  // in Gauss
-
-  for (int i = 0; i<MetaData.TopGridRank; i++) {ConductionCloudVelocity[i] = 0.0;}
- 
-  float ConductionCloudPulseHeight;
-  FLOAT ConductionCloudPulseWidth;
-  int ConductionCloudPulseType = 0;
-
-  TestProblemData.MultiSpecies = MultiSpecies;  // set this from global data (kind of a hack, but necessary)
-
-  // Read parameters
-  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
-    ret = 0;
-    ret += sscanf(line, "ConductionCloudPulseHeight = %"FSYM, &ConductionCloudPulseHeight);
-    ret += sscanf(line, "ConductionCloudTotalEnergy = %"FSYM, &ConductionCloudTotalEnergy);
-    ret += sscanf(line, "ConductionCloudTemperature = %"FSYM, &ConductionCloudTemperature);
-    ret += sscanf(line, "ConductionCloudDensity = %"FSYM, &ConductionCloudDensity);
-    ret += sscanf(line, "ConductionCloudPulseWidth = %"PSYM, &ConductionCloudPulseWidth);
-    ret += sscanf(line, "ConductionCloudPulseType = %"ISYM, &ConductionCloudPulseType);
-
-    ret += sscanf(line, "TestProblemUseMetallicityField  = %"ISYM, &TestProblemData.UseMetallicityField);
-    ret += sscanf(line, "TestProblemInitialMetallicityFraction  = %"FSYM, &TestProblemData.MetallicityField_Fraction);
-
-    /* read in more general test parameters to set species, turn on color fields, etc. */
-    ret += sscanf(line, "TestProblemHydrogenFractionByMass = %"FSYM, &TestProblemData.HydrogenFractionByMass);
-    ret += sscanf(line, "TestProblemDeuteriumToHydrogenRatio = %"FSYM, &TestProblemData.DeuteriumToHydrogenRatio);
-
-    ret += sscanf(line, "TestProblemInitialHIFractionInner  = %"FSYM, &TestProblemData.HI_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHIIFractionInner  = %"FSYM, &TestProblemData.HII_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHeIFractionInner  = %"FSYM, &TestProblemData.HeI_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHeIIFractionInner  = %"FSYM, &TestProblemData.HeII_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHeIIIFractionInner  = %"FSYM, &TestProblemData.HeIII_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHMFractionInner  = %"FSYM, &TestProblemData.HM_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialH2IFractionInner  = %"FSYM, &TestProblemData.H2I_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialH2IIFractionInner  = %"FSYM, &TestProblemData.H2II_Fraction_Inner);
-
-    ret += sscanf(line, "TestProblemInitialDIFractionInner  = %"FSYM, &TestProblemData.DI_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialDIIFractionInner  = %"FSYM, &TestProblemData.DII_Fraction_Inner);
-    ret += sscanf(line, "TestProblemInitialHDIFractionInner  = %"FSYM, &TestProblemData.HDI_Fraction_Inner);
-
-    ret += sscanf(line, "TestProblemInitialHIFraction  = %"FSYM, &TestProblemData.HI_Fraction);
-    ret += sscanf(line, "TestProblemInitialHIIFraction  = %"FSYM, &TestProblemData.HII_Fraction);
-    ret += sscanf(line, "TestProblemInitialHeIFraction  = %"FSYM, &TestProblemData.HeI_Fraction);
-    ret += sscanf(line, "TestProblemInitialHeIIFraction  = %"FSYM, &TestProblemData.HeII_Fraction);
-    ret += sscanf(line, "TestProblemInitialHeIIIFraction  = %"FSYM, &TestProblemData.HeIII_Fraction);
-    ret += sscanf(line, "TestProblemInitialHMFraction  = %"FSYM, &TestProblemData.HM_Fraction);
-    ret += sscanf(line, "TestProblemInitialH2IFraction  = %"FSYM, &TestProblemData.H2I_Fraction);
-    ret += sscanf(line, "TestProblemInitialH2IIFraction  = %"FSYM, &TestProblemData.H2II_Fraction);
-
-    ret += sscanf(line, "TestProblemInitialDIFraction  = %"FSYM, &TestProblemData.DI_Fraction);
-    ret += sscanf(line, "TestProblemInitialDIIFraction  = %"FSYM, &TestProblemData.DII_Fraction);
-    ret += sscanf(line, "TestProblemInitialHDIFraction  = %"FSYM, &TestProblemData.HDI_Fraction);
-
-    ret += sscanf(line, "TestProblemUseMetallicityField  = %"ISYM, &TestProblemData.UseMetallicityField);
-    ret += sscanf(line, "TestProblemInitialMetallicityFraction  = %"FSYM, &TestProblemData.MetallicityField_Fraction);
-
-    ret += sscanf(line, "TestProblemMultiMetals  = %"ISYM, &TestProblemData.MultiMetals);
-    ret += sscanf(line, "TestProblemInitialMultiMetalsField1Fraction  = %"FSYM, &TestProblemData.MultiMetalsField1_Fraction);
-    ret += sscanf(line, "TestProblemInitialMultiMetalsField2Fraction  = %"FSYM, &TestProblemData.MultiMetalsField2_Fraction);
-
-    if (ret == 0 && 
-	strstr(line, "=") && strstr(line, "ConductionCloud") &&
-	line[0] != '#' && MyProcessorNumber == ROOT_PROCESSOR) {
-      fprintf(stderr, "*** warning: the following parameter line was not interpreted:\n%s\n", line);
-    }
-  }
-
-  float DensityUnits=1.0, LengthUnits=1.0, TemperatureUnits=1.0, TimeUnits=1.0,
-    VelocityUnits=1.0;
-  double MassUnits=1.0;
-
-  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
-	       &TimeUnits, &VelocityUnits, &MassUnits, 0.0) == FAIL) {
-    fprintf(stderr, "Error in GetUnits.\n");
-    return FAIL;
-  }
-
-  float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
-
-  if(ConductionCloudTemperature > 1.0){
-
-    ConductionCloudTotalEnergy = (Boltzmann*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
-    ConductionCloudTotalEnergy /= (VelocityUnits*VelocityUnits);
-    printf("ConductionCloudTotalEnergy is %e and ConductionCloudTemperature is %e\n\n",ConductionCloudTotalEnergy, ConductionCloudTemperature);
-    fflush(stdout);
-  }
-
-  // Create a uniform grid
-  if (TopGrid.GridData->InitializeUniformGrid(ConductionCloudDensity,
-					      ConductionCloudTotalEnergy,
-					      ConductionCloudTotalEnergy,
-					      ConductionCloudVelocity,
-					      ConductionCloudInitialUniformBField) == FAIL) {
-    ENZO_FAIL("Error in InitializeUniformGrid.");
-  }
-  
-  // Then perturb it
-  if (TopGrid.GridData->ConductionCloudInitialize(ConductionCloudPulseHeight, 
-						 ConductionCloudPulseWidth, 
-						 ConductionCloudPulseType) == FAIL) {
-    ENZO_FAIL("Error in ConductionCloudInitialize.");
-  }
-
-  // set up field names and units
-  i = 0;
-  DataLabel[i++] = "Density";
-  DataLabel[i++] = "TotalEnergy";
-  if (DualEnergyFormalism) {DataLabel[i++] = "GasEnergy";}
-  if (MetaData.TopGridRank > 0) {DataLabel[i++] = "x-velocity";}
-  if (MetaData.TopGridRank > 1) {DataLabel[i++] = "y-velocity";}
-  if (MetaData.TopGridRank > 2) {DataLabel[i++] = "z-velocity";}
-
-  if (TestProblemData.MultiSpecies) {
-    DataLabel[i++] = "Electron_Density";
-    DataLabel[i++] = "HI_Density";
-    DataLabel[i++] = "HII_Density";
-    DataLabel[i++] = "HeI_Density";
-    DataLabel[i++] = "HeII_Density";
-    DataLabel[i++] = "HeIII_Density";
-    if (TestProblemData.MultiSpecies > 1) {
-      DataLabel[i++] = "HM_Density";
-      DataLabel[i++] = "H2I_Density";
-      DataLabel[i++] = "H2II_Density";
-    }
-    if (TestProblemData.MultiSpecies > 2) {
-      DataLabel[i++] = "DI_Density";
-      DataLabel[i++] = "DII_Density";
-      DataLabel[i++] = "HDI_Density";
-    }
-  }
-
-  if (TestProblemData.UseMetallicityField)
-    DataLabel[i++] = "Metal_Density";
-
-  for (j=0; j < i; j++) 
-    DataUnits[j] = NULL;
-
-  /* Write parameters to parameter output file */
- 
-  if (MyProcessorNumber == ROOT_PROCESSOR) {
-    fprintf(Outfptr, "ConductionCloudPulseHeight = %"FSYM"\n", ConductionCloudPulseHeight);
-    fprintf(Outfptr, "ConductionCloudPulseWidth = %"PSYM"\n", ConductionCloudPulseWidth);
-    fprintf(Outfptr, "ConductionCloudPulseType = %"ISYM"\n", ConductionCloudPulseType);
-
-    fprintf(Outfptr, "TestProblemUseMetallicityField  = %"ISYM"\n", TestProblemData.UseMetallicityField);
-    fprintf(Outfptr, "TestProblemInitialMetallicityFraction  = %"FSYM"\n", TestProblemData.MetallicityField_Fraction);
-
-
-    fprintf(Outfptr, "TestProblemHydrogenFractionByMass = %"FSYM"\n",   TestProblemData.HydrogenFractionByMass);
-    fprintf(Outfptr, "TestProblemDeuteriumToHydrogenRatio = %"FSYM"\n", TestProblemData.DeuteriumToHydrogenRatio);
-
-    fprintf(Outfptr, "TestProblemInitialHIFractionInner  = %"FSYM"\n", TestProblemData.HI_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHIIFractionInner  = %"FSYM"\n", TestProblemData.HII_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHeIFractionInner  = %"FSYM"\n", TestProblemData.HeI_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHeIIFractionInner  = %"FSYM"\n", TestProblemData.HeII_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHeIIIFractionInner  = %"FSYM"\n", TestProblemData.HeIII_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHMFractionInner  = %"FSYM"\n", TestProblemData.HM_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialH2IFractionInner  = %"FSYM"\n", TestProblemData.H2I_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialH2IIFractionInner  = %"FSYM"\n", TestProblemData.H2II_Fraction_Inner);
-
-    fprintf(Outfptr, "TestProblemInitialDIFractionInner  = %"FSYM"\n", TestProblemData.DI_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialDIIFractionInner  = %"FSYM"\n", TestProblemData.DII_Fraction_Inner);
-    fprintf(Outfptr, "TestProblemInitialHDIFractionInner  = %"FSYM"\n", TestProblemData.HDI_Fraction_Inner);
-
-    fprintf(Outfptr, "TestProblemInitialHIFraction  = %"FSYM"\n", TestProblemData.HI_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHIIFraction  = %"FSYM"\n", TestProblemData.HII_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHeIFraction  = %"FSYM"\n", TestProblemData.HeI_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHeIIFraction  = %"FSYM"\n", TestProblemData.HeII_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHeIIIFraction  = %"FSYM"\n", TestProblemData.HeIII_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHMFraction  = %"FSYM"\n", TestProblemData.HM_Fraction);
-    fprintf(Outfptr, "TestProblemInitialH2IFraction  = %"FSYM"\n", TestProblemData.H2I_Fraction);
-    fprintf(Outfptr, "TestProblemInitialH2IIFraction  = %"FSYM"\n", TestProblemData.H2II_Fraction);
-
-    fprintf(Outfptr, "TestProblemInitialDIFraction  = %"FSYM"\n", TestProblemData.DI_Fraction);
-    fprintf(Outfptr, "TestProblemInitialDIIFraction  = %"FSYM"\n", TestProblemData.DII_Fraction);
-    fprintf(Outfptr, "TestProblemInitialHDIFraction  = %"FSYM"\n", TestProblemData.HDI_Fraction);
-
-    fprintf(Outfptr, "TestProblemUseMetallicityField  = %"ISYM"\n", TestProblemData.UseMetallicityField);
-    fprintf(Outfptr, "TestProblemInitialMetallicityFraction  = %"FSYM"\n", TestProblemData.MetallicityField_Fraction);
-