Commits

Brian O'Shea  committed 3324985 Merge

merged from tip.

  • Participants
  • Parent commits 0ea9dba, 1332329
  • Branches rotating-sphere

Comments (0)

Files changed (53)

File doc/manual/source/EnzoParameters.rst

     Causes the simulation to immediately stop when a specified level is
     reached. Default value 0 (off), possible values are levels 1
     through maximum number of levels in a given simulation.
+``NumberOfOutputsBeforeExit`` (external)
+    After this many datadumps have been written, the code will exit.  If 
+    set to 0 (default), this option will not be used.  Default: 0.
 ``StopCPUTime`` (external)
     Causes the simulation to stop if the wall time exceeds ``StopCPUTime``.
     The simulation will output if the wall time after the next
     is valid for temperatures greater than 10,000 K. This requires the
     file ``TREECOOL`` to execute. Default: 0
 ``MetalCooling`` (external)
-    This flag (0 - off, 1 - metal cooling from Glover & Jappsen 2007, 2
-    - Cen, 3 - Cloudy cooling from Smith, Sigurdsson, & Abel 2008)
-    turns on metal cooling for runs that track metallicity. Option 1 is
-    valid for temperatures between 100 K and 10\ :sup:`8`\  K because
-    it considers fine-structure line emission from carbon, oxygen, and
-    silicon and includes the additional metal cooling rates from
-    Sutherland & Dopita (1993). Option 2 is only valid for temperatures
-    above 10\ :sup:`4`\  K. Option 3 uses multi-dimensional tables of
-    heating/cooling values created with Cloudy and optionally coupled
-    to the ``MultiSpecies`` chemistry/cooling solver. This method is valid
-    from 10 K to 10\ :sup:`8`\  K. See the Cloudy Cooling parameters below.
-    Default: 0.
+    This flag (0 - off, 1 - metal cooling from Glover & Jappsen 2007,
+    2 - Cen et al (1995), 3 - Cloudy cooling from Smith, Sigurdsson, &
+    Abel 2008) turns on metal cooling for runs that track
+    metallicity. Option 1 is valid for temperatures between 100 K and
+    10\ :sup:`8`\ K because it considers fine-structure line emission
+    from carbon, oxygen, and silicon and includes the additional metal
+    cooling rates from Sutherland & Dopita (1993). Option 2 is only
+    valid for temperatures above 10\ :sup:`4`\ K. Option 3 uses
+    multi-dimensional tables of heating/cooling values created with
+    Cloudy and optionally coupled to the ``MultiSpecies``
+    chemistry/cooling solver. This method is valid from 10 K to 10\
+    :sup:`8`\ K. See the Cloudy Cooling parameters below.  Default: 0.
 ``MetalCoolingTable`` (internal)
     This field contains the metal cooling table required for
     ``MetalCooling`` option 1. In the top level directory input/, there are
     hydrogen (H2) dissociation rate. There a normalization is performed
     on the rate by multiplying it with ``RadiationSpectrumNormalization``.
     Default: 1e-21
+``RadiationFieldRedshift`` (external)
+    This parameter specifies the redshift at which the radiation field
+    is calculated.  Default: 0
 ``RadiationShield`` (external)
     This parameter specifies whether the user wants to employ
     approximate radiative-shielding. This parameter will be
 
 ``UseMinimumPressureSupport`` (external)
     When radiative cooling is turned on, and objects are allowed to
-    collapse to very small sizes (i.e. a few cells), and they are
-    evolved for many, many dynamical times, then unfortunate things
-    happen. Primarily, there is some spurious angular momentum
-    generation, and possible some resulting momentum non-conservation.
-    To alleviate this problem, a very simple fudge was introduced: if
-    this flag is turned on, then a minimum temperature is applied to
-    grids with level == ``MaximumRefinementLevel``. This minimum
-    temperature is that required to make each cell Jeans stable
-    multiplied by the parameter below. If you use this, it is advisable
-    to also set the gravitational smoothing length in the form of
-    ``MaximumGravityRefineLevel`` to 2 or 3 less than
-    ``MaximumRefinementLevel``. Default: 0
+    collapse to very small sizes so that their Jeans length is no
+    longer resolved, then they may undergo artificial fragmentation
+    and angular momentum non-conservation.  To alleviate this problem,
+    as discussed in more detail in Machacek, Bryan & Abel (2001), a
+    very simple fudge was introduced: if this flag is turned on, then
+    a minimum temperature is applied to grids with level ==
+    ``MaximumRefinementLevel``. This minimum temperature is that
+    required to make each cell Jeans stable multiplied by the
+    parameter below.  More precisely, the temperature of a cell is set
+    such that the resulting Jeans length is the square-root of the
+    parameter ``MinimumPressureSupportParameter``.  So, for the
+    default value of 100 (see below), this insures that the ratio of
+    the Jeans length/cell size is at least 10.  Default: 0
 ``MinimumPressureSupportParameter`` (external)
-    This is the parameter alluded to above. Very roughly speaking, is
-    the number of cells over which a gravitationally bound small cold
-    clump, on the most refined level, will be spread over. Default:
-    100
+    This is the numerical parameter discussed above. Default: 100
 
 Radiative Transfer (Ray Tracing) Parameters
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     ``mbh_particle_io.dat``
 
 Conduction
-^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~~~~~~~
 
 Isotropic and anisotropic thermal conduction are implemented using the
 method of Parrish and Stone: namely, using an explicit, forward
     a value of 0.5 or less.
     Default: 0.5
 
+Shock Finding Parameters
+~~~~~~~~~~~~~~~~~~~~~~~~
+For details on shock finding in Enzo see :ref:`shock_finding`.
+
+``ShockMethod`` (external)
+    This parameter controls the use and type of shock finding. Default: 0
+    
+    ::
+
+	0 - Off
+	1 - Temperature Dimensionally Unsplit Jumps
+	2 - Temperature Dimensionally Split Jumps
+	1 - Velocity Dimensionally Unsplit Jumps
+	2 - Velocity Dimensionally Split Jumps
+
+``ShockTemperatureFloor`` (external)
+    When calculating the mach number using temperature jumps, set the
+    temperature floor in the calculation to this value.
+
+``StorePreShockFields`` (external)
+    Optionally store the Pre-shock Density and Temperature during data output.
+
+
 .. _testproblem_param:
 
 Test Problem Parameters

File doc/manual/source/physics/cooling.rst

 
    ``MetalCooling`` = 1
 
-#. Cen cooling.
+#. Cen et al (1995) cooling. This uses output from a Raymond-Smith
+   code to determine cooling rates from T > 10\ :sup:`4`\ K.  No ionizing
+   background is used in computing cooling rates.  This method has
+   not been extensively tested in the context of Enzo.
 
    ``RadiativeCooling`` = 1
 
 dependent and work by setting the photoionization and photoheating
 coeffiecients for H, He, and He\ :sup:`+`\.  See
 :ref:`radiation_backgrounds` for the additional parameters that
-control the UV backgrounds.
+control the UV backgrounds.

File doc/manual/source/physics/cosmic_rays.rst

-.. _cosmic_rays:
-
-Cosmic Rays
-==================
-.. sectionauthor:: Sam Skillman <Samuel.Skillman@Colorado.EDU>
-.. versionadded:: 2.1
-

File doc/manual/source/physics/index.rst

    hydro_methods.rst
    cooling.rst
    radiative_transfer.rst
-   cosmic_rays.rst
    shock_finding.rst
 

File doc/manual/source/physics/shock_finding.rst

 
 Shock Finding
 ==================
-.. sectionauthor:: Sam Skillman <Samuel.Skillman@Colorado.EDU>
+.. sectionauthor:: Sam Skillman <samskillman@gmail.com>
 .. versionadded:: 2.1
+*Source: Grid_FindShocks.C*
 
+Shock finding is accomplished using one of four methods.  The primary
+method uses a coordinate-unsplit temperature jump (method 1), as described in
+`Skillman et. al. 2008
+<http://adsabs.harvard.edu/abs/2008ApJ...689.1063S>`_ with the
+exception that instead of searching across multiple grids for the pre-
+and post-shock cells, we terminate the search at the edge of the ghost
+zones within each grid.  
+
+Shock finding is controlled by the ``ShockMethod`` parameter, which
+can take the following values:
+
+0 - Off
+
+1 - Unsplit Temperature Jumps
+
+2 - Dimensionally Split Temperature Jumps
+
+3 - Unsplit Velocity Jumps
+
+4 - Dimensionally Split Velocity Jumps
+
+When ``ShockMethod`` nonzero, this will create a "Mach" field in the
+output files.   
+
+Note: Method 1 has been used the most by the developer, and therefore
+is the primary method.  Method 2 has been tested quite a bit, but the
+downsides of using a dimensionally split method are outlined in the
+above paper.  Methods 3 and 4 are more experimental and will run, but
+results may vary.
+
+Additional Shock Finding Parameters:
+
+``ShockTemperatureFloor`` - When calculating the mach number using temperature jumps, set the temperature floor in the calculation to this value.
+
+``StorePreShockFields`` - Optionally store the Pre-shock Density and Temperature during data output.
+
+
+
+
+

File doc/manual/source/physics/star_particles.rst

      respectively. Note that :math:`\rho_{\mathrm{th}}` is not the same as the
      critical density for star formation listed above. :math:`\rho` is the
      gas density of the cell.
+
      .. math::
-     
-        t_{\ast}(\rho)=t_0^{\ast}\left(\frac{\rho}{\rho_{\mathrm{th}}}\right)^{-1/2}
+
+     t_{\ast}(\rho)=t_0^{\ast}\left(\frac{\rho}{\rho_{\mathrm{th}}}\right)^{-1/2}
   
 * Mass fraction in cold clouds, :math:`x` (see Eqns. 16 and 18).
      :math:`y` is a dimensionless quantity
      :math:`\epsilon_{\textrm{SN}}` is the energy released from a nominal
      supernova and is set to 4e48 ergs; and finally :math:`\Lambda(\rho, T, z)`
      is the cooling rate of the cell of gas.
+
      .. math::
      
         y\equiv\frac{t_{\ast}\Lambda(\rho,T,z)}{\rho[\beta u_{\mathrm{SN}}-(1-\beta)u_{\mathrm{SN}}]}
 :math:`m` is the mass of gas in the cell;
 :math:`\Delta t` is the size of the simulation time step that
 is operative for the cell (which changes over AMR levels, of course).
+
 .. math::
 
-   p_{\ast}=\frac{m}{m_{\ast}}\left\{1-\exp\left[-\frac{(1-\beta)x\Delta t}{t_{\ast}}\right]\right\}\
+   p_{\ast}=\frac{m}{m_{\ast}}\left\{1-\exp\left[-\frac{(1-\beta)x\Delta t}{t_{\ast}}\right]\right\}
 
 If this star formula is used with AMR, some caution is required. Primarily,
 the AMR refinement can not be too aggressive. Values of ``OverDensityThreshold``

File doc/manual/source/user_guide/RunningLargeSimulations.rst

 
 Here we describe how to efficiently run a large simulation on a high
 number of processors, such as particular parameters to set and
-suggested number of MPI tasks for a given problem size.
+suggested number of MPI tasks for a given problem size.  For a problem
+to be scalable, most of the code must be parallel to achieve high
+performance numbers on large MPI process counts (see `Amdahl's
+Law`__).  In general, the user wants to pick the number of processors
+so that computation is still dominant over communication time.  If the
+processor count is too high, communication time will become too large
+and might even *slow* down the simulation!
+
+For picking the number of processors for an Enzo run, a good starting
+point is putting a 64\ :sup:`3` box on each processor for both AMR and
+unigrid setups.  For example, a 256\ :sup:`3` simulation would run
+well on (256/64)\ :sup:`3` = 64 processors.  For nested grid
+simulations, the outer boxes usually require little computation
+compared to the "zoom-in" region, so the processor count should be
+based on the inner-most nested grid size.  The user can experiment
+with increasing the processor count from this suggestion, but strong
+scaling (i.e. linear speedup with processor count) is not to be
+expected.  Little performance gains (as of v2.0) can be expected
+beyond assigning a 32\ :sup:`3` cube per processor.
+
+.. note:: 
+
+   The level-0 grid is only partitioned during the problem
+   initialization.  It will *never* be re-partitioned if the user
+   restarts with a different number of processors.  However, some
+   performance gains can be expected even if a processor does not
+   contain a level-0 grid because of the work on finer levels.
+
+Important Parameters
+--------------------
+
+* ``LoadBalancing``: Default is 1, which moves work from overloaded to
+  underutilized processes, regardless of the grid position.  **New for
+  v2.1**: In some cases but not always, speedups can be found in load
+  balancing on a `space filling curve`_ (``LoadBalancing = 4``).  Here
+  the grids on each processor will be continuous on the space filling
+  curve.  This results in a grouped set of grids, requiring less
+  communication from other processors (and even other compute nodes).
+
+* ``SubgridSizeAutoAdjust`` and ``OptimalSubgridsPerProcessor``: **New for
+  v2.1** Default is ON and 16, respectively.  The maximum subgrid size
+  and edge length will be dynamically adjusted on each AMR level
+  according to the number of cells on the level and number of
+  processors.  The basic idea behind increasing the subgrid sizes
+  (i.e. coalescing grids) is to reduce communication between grids.
+
+* ``MinimumSubgridEdge`` and ``MaximumSubgridSize``: *Unused if
+  SubgridAutoAdjust is ON*.  Increase both of these parameters to
+  increase the average subgrid size, which might reduce communication
+  and speedup the simulation.
+
+* ``UnigridTranspose``: Default is 0, which is employs blocking MPI
+  communication to transpose the root grid before and after the FFT.
+  In level-0 grids |ge| 1024\ :sup:`3`, this becomes the most
+  expense part of the calculation.  In these types of large runs,
+  Option 2 is recommended, which uses non-blocking MPI calls; however
+  it has some additional memory overhead, which is the reason it is
+  not used by default.
+
+Compile-time options
+--------------------
+
+* ``max-subgrids``: If the number of subgrids in a single AMR level
+  exceeds this value, then the simulation will crash.  Increase as
+  necessary.  Default: 100,000
+
+* ``ooc-boundary-yes``: Stores the boundary conditions out of core,
+  i.e. on disk.  Otherwise, each processor contains a complete copy of
+  the external boundary conditions.  This becomes useful in runs with
+  large level-0 grids.  For instance in a 1024\ :sup:`3` simulation
+  with 16 baryon fields, each processor will contain a set of
+  boundary conditions on 6 faces of 1024\ :sup:`2` with 16 baryon
+  fields.  In single precision, this requires 402MB!  Default: OFF
+
+* ``fastsib-yes``: Uses a chaining mesh to help locate sibling grids
+  when constructing the boundary conditions.  Default: ON
+
+.. |ge| unicode:: 0x2265
+
+.. _space filling curve: http://en.wikipedia.org/wiki/Hilbert_curve
+
+.. __: http://en.wikipedia.org/wiki/Amdahl's_law

File doc/manual/source/user_guide/index.rst

    EmbeddedPython.rst
    HierarchyFile.rst
    FlowChart.rst
-   Presentations.rst

File run/test_makespreadsheet.py

+#!/usr/bin/env python
+#
+#  Looks through the run directory for "*.enzotest" and makes a csv spreadsheet of all answer test scripts and their properties.
+#
+#  Author: David Collins (dcollins4096@gmail.com), 2011-06-14 11:19 AM.  It's a bright sunny day here in Los Alamos.
+#
+
+import fnmatch
+import os
+
+
+#Hunt for enzotest files.
+dbg = 0
+matches = []
+for root, dirnames, filenames in os.walk('.'):
+  for filename in fnmatch.filter(filenames, '*.enzotest'):
+      matches.append(os.path.join(root, filename))
+
+#Generate dictionary, make list of attributes
+tests={}
+attribute_list=['name','nprocs','max_time_minutes','dimensionality','runtime','critical','cadence','answer_testing_script','hydro','gravity','cooling','chemistry','cosmology','author','mhd','radiation','AMR']
+for file in matches:
+    if dbg > 0:
+        print file
+    lines = open(file,'r').readlines()
+    tests[file]={}
+    for line in lines:
+        if line.strip():
+            key, value = line.split("=")
+            if value[-1]=='\n':
+                value = value[:-1]
+            if value.count("#") > 0:
+                value = value[0:value.index("#")]
+            key=key.strip()
+            value=value.strip()
+            if key not in attribute_list:  
+                attribute_list.append(key)
+            tests[file][key]=value
+
+#make csv
+csv = open('test_spreadsheet.csv','w')
+dir( csv)
+
+#head row
+csv.write("filename"+",%s"*len(attribute_list)%tuple(attribute_list)+"\n")
+
+#lines
+for file in matches:
+    csv.write(file)
+    for attr in attribute_list:
+        csv.write(", %s"%(tests[file].get(attr,'')))
+    csv.write("\n")
+csv.close()
+
+
+#end

File src/enzo/CheckForOutput.C

   }
 #endif   
 
+  if (MetaData.NumberOfOutputsBeforeExit && MetaData.WroteData) {
+    MetaData.OutputsLeftBeforeExit--;
+    if (MetaData.OutputsLeftBeforeExit <= 0) {
+      if (MyProcessorNumber == ROOT_PROCESSOR) {
+	fprintf(stderr, "Exiting after writing %"ISYM" datadumps.\n",
+		MetaData.NumberOfOutputsBeforeExit);
+      }      
+      my_exit(EXIT_SUCCESS);
+    }
+  }
+
   return SUCCESS;
 }

File src/enzo/CommunicationUpdateStarParticleCount.C

  
   delete [] TotalParticleCount;
   delete [] PartialParticleCount;
+  delete [] TotalStarParticleCount;
+  delete [] PartialStarParticleCount;
  
   LCAPERF_STOP("UpdateStarParticleCount");
  
  
   delete [] TotalParticleCount;
   delete [] PartialParticleCount;
+  delete [] TotalStarParticleCount;
+  delete [] PartialStarParticleCount;
 
   return SUCCESS;
 }

File src/enzo/CosmicRayData.h

-/***********************************************************************
-/
-/  Cosmic Ray Acceleration Efficiency Data		
-/
-***********************************************************************/
-
-struct CosmicRayDataType
-{
-  // Number of dimension to interpolate over
-  int CosmicRayGridRank;
-
-  // Maximum number of values for Pre-CR population
-  int CRMaxNumberPrePopValues;
-  // Maximum number of values for Mach 
-  int CRMaxNumberMachValues;
-
-  // Cooling grid run file
-  char *CosmicRayGridRunFile;
-
-  // Values for first parameter space (Pre-Shock CR Percentage)
-  float *CRPrePopValues;
-  // Values for second parameter space (Mach Number)
-  float *CRMachValues;
-
-  //Min/max of the paramters
-  float CRMinPrePop;
-  float CRMaxPrePop;
-  float CRMinMach;
-  float CRMaxMach;
-
-  // Number of values for first parameter space
-  int CRNumberPrePopValues;
-  // Number of values for second parameter space
-  int CRNumberMachValues;
-
-  // Efficiency Values
-  float **CREfficiency;
-};

File src/enzo/CosmologySimulationInitialize.C

   char *GPotName  = "Grav_Potential";
   char *ForbidName  = "ForbiddenRefinement";
   char *MachName   = "Mach";
-  char *CRName     = "CR_Density";
   char *PSTempName = "PreShock_Temperature";
   char *PSDenName  = "PreShock_Density";
   char *ExtraNames[2] = {"Z_Field1", "Z_Field2"};
   if (WritePotential)
     DataLabel[i++] = GPotName;  
 
-  if (CRModel) {
+#ifdef EMISSIVITY
+  if (StarMakerEmissivityField > 0)
+    DataLabel[i++] = EtaName;
+#endif
+ 
+  if (ShockMethod) {
     DataLabel[i++] = MachName;
     if(StorePreShockFields){
       DataLabel[i++] = PSTempName;
       DataLabel[i++] = PSDenName;
     }
-    DataLabel[i++] = CRName;
   } 
 
-#ifdef EMISSIVITY
-  if (StarMakerEmissivityField > 0)
-    DataLabel[i++] = EtaName;
-#endif
- 
   for (j = 0; j < i; j++)
     DataUnits[j] = NULL;
  

File src/enzo/EvolveLevel.C

 /                Cleaned up error handling and created new routines for
 /                computing the timestep, output, handling fluxes
 /  modified10: July, 2009 by Sam Skillman
-/                Added shock and cosmic ray analysis
+/                Added shock analysis
 /
 /  PURPOSE:
 /    This routine is the main grid evolution function.  It assumes that the
       Grids[grid1]->GridData->StarParticleHandler
 	(Grids[grid1]->NextGridNextLevel, level ,dtLevelAbove);
 
-      /* Include shock-finding and cosmic ray acceleration */
+      /* Include shock-finding */
 
       Grids[grid1]->GridData->ShocksHandler();
 

File src/enzo/EvolvePhotons.C

 	  else
 	    Helper = NULL;
 
+#ifdef BITWISE_IDENTICALITY
+	  Temp->GridData->PhotonSortLinkedLists();
+#endif
 	  Temp->GridData->TransportPhotonPackages
 	    (lvl, &PhotonsToMove, GridNum, Grids0, nGrids0, Helper, 
 	     Temp->GridData);

File src/enzo/Grid.h

   FLOAT *CellLeftEdge[MAX_DIMENSION];
   FLOAT *CellWidth[MAX_DIMENSION];
   fluxes *BoundaryFluxes;
+  float *YT_TemperatureField;                         // place to store temperature field
+                                                      // for call to yt.
 
   // For restart dumps
 
    float GadgetCoolingRateFromU(float u, float rho, float *ne_guess, 
 				float redshift);
 
-// Functions for shock finding and cosmic ray acceleration
+// Functions for shock finding
 //
    int FindShocks();
    int FindTempSplitShocks();
 			    int &HMNum, int &H2INum, int &H2IINum,
                             int &DINum, int &DIINum, int &HDINum);
 
-/* Identify shock/cosmic ray fields. */
-  int IdentifyCRSpeciesFields(int &MachNum, int&CRNum, 
-			      int &PSTempNum, int &PSDenNum);
+  /* Identify shock fields. */
+  int IdentifyShockSpeciesFields(int &MachNum,int &PSTempNum, int &PSDenNum);
+
   // Identify Simon Glover Species Fields
   int IdentifyGloverSpeciesFields(int &HIINum,int &HINum,int &H2INum,
 				  int &DINum,int &DIINum,int &HDINum,

File src/enzo/Grid_ComputeAccelerationFieldExternal.C

 
 	  
 	  radius = sqrt(rsquared);
-	  rcore = PointSourceGravityCoreRadius/LengthUnits;
+	  
+	  if(ProblemType == 31){
+	    rcore = PointSourceGravityCoreRadius;  // already in code units                                                      
+	  } else {
+	    rcore = PointSourceGravityCoreRadius/LengthUnits;  // convert from CGS to code                                       
+	  }
 	  FLOAT x = radius/rcore;
-	  accel = GravConst*PointSourceGravityConstant*
+	  accel = GravConst*PointSourceGravityConstant*MassUnitsDouble*
 	    ((log(1.0+x  )-x  /(1.0+x  )) /
 	     (log(1.0+1.0)-1.0/(1.0+1.0))) / 
 	    POW(radius*LengthUnits, 2) / AccelUnits;

File src/enzo/Grid_ComputePressure.C

 	BaryonField[TENum][i] = e + kineticE + OneHalf*B2/density;
 
       } else { 
-      /* gas energy = E - 1/2 v^2. */
-	gas_energy    = total_energy - OneHalf*(velocity1*velocity1 +
-						velocity2*velocity2 +
-						velocity3*velocity3);
+	if (DualEnergyFormalism == 0){ 
+	  /* gas energy = E - 1/2 v^2. */
+	  gas_energy    = total_energy - OneHalf*(velocity1*velocity1 +
+						  velocity2*velocity2 +
+						  velocity3*velocity3);
+	} else {
+	  gas_energy = BaryonField[GENum][i];
+	}
+
  	if (HydroMethod == MHD_RK) {
 	  float B2 = pow(BaryonField[B1Num][i],2) + pow(BaryonField[B2Num][i],2) +
 	    pow(BaryonField[B3Num][i],2);

File src/enzo/Grid_ConvertToNumpy.C

 
 	/* Get grid temperature field. */
 	int size = 1;
-	float *temperature;
 	for (dim = 0; dim < GridRank; dim++)
 	  size *= GridDimension[dim];
-	temperature = new float[size];
-	if (this->ComputeTemperatureField(temperature) == FAIL) {
+	if (YT_TemperatureField != NULL) {
+	  delete [] YT_TemperatureField;
+	  YT_TemperatureField = NULL;
+	}
+	YT_TemperatureField = new float[size];
+	if (this->ComputeTemperatureField(YT_TemperatureField) == FAIL) {
 	  ENZO_FAIL("Error in grid->ComputeTemperatureField.\n");
 	}
 	dataset = (PyArrayObject *) PyArray_SimpleNewFromData(
-	        3, dims, ENPY_BFLOAT, temperature);
+	        3, dims, ENPY_BFLOAT, YT_TemperatureField);
 	dataset->flags &= NPY_OWNDATA;
 	PyDict_SetItemString(grid_data, "Temperature", (PyObject*) dataset);
 	Py_DECREF(dataset);

File src/enzo/Grid_CopyParentToGravitatingFieldBoundary.C

      the region dim (in parent units). */
  
   for (dim = 0; dim < GridRank; dim++) {
-    // SubGridExtra[dim] = nint((GridLeftEdge[dim] -
-    // 			      GravitatingMassFieldLeftEdge[dim])
-    // 			     /GravitatingMassFieldCellSize);
-       SubGridExtra[dim] = nint((CellLeftEdge[dim][0] -
-    			      GravitatingMassFieldLeftEdge[dim])
+    SubGridExtra[dim] = nint((GridLeftEdge[dim] -
+			      GravitatingMassFieldLeftEdge[dim])
     			     /GravitatingMassFieldCellSize);
+    // SubGridExtra[dim] = nint((CellLeftEdge[dim][0] -
+    //  			      GravitatingMassFieldLeftEdge[dim])
+    //  			     /GravitatingMassFieldCellSize);
     ParentOffset[dim] = nint((GravitatingMassFieldLeftEdge[dim] -
 		  ParentGrid->GravitatingMassFieldLeftEdge[dim])/
 			     GravitatingMassFieldCellSize);

File src/enzo/Grid_CosmologySimulationInitializeGrid.C

 #ifdef EMISSIVITY
   int EtaNum;
 #endif
-  int CRNum, MachNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
  
   int ExtraField[2];
   int ForbidNum, iTE;
     if (StarMakerEmissivityField > 0)
       FieldType[EtaNum = NumberOfBaryonFields++] = Emissivity0;
 #endif
-    if(CRModel){
+    if(ShockMethod){
       FieldType[MachNum   = NumberOfBaryonFields++] = Mach;
       if(StorePreShockFields){
 	FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature;
 	FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity;
       }
-      FieldType[CRNum     = NumberOfBaryonFields++] = CRDensity;
     }    
   }
  
       }
  
       //Shock/Cosmic Ray Model
-      if (CRModel && ReadData) {
+      if (ShockMethod && ReadData) {
 	BaryonField[MachNum][i] = tiny_number;
-	BaryonField[CRNum][i] = tiny_number;
 	if (StorePreShockFields) {
 	  BaryonField[PSTempNum][i] = tiny_number;
 	  BaryonField[PSDenNum][i] = tiny_number;

File src/enzo/Grid_FindShocks.C

 /***********************************************************************
 /
-/  GRID CLASS (Find Shocks and Accelerate Cosmic Rays)
+/  GRID CLASS (Find Shocks)
 /
 /  written by: Sam Skillman
 /  date:       May, 2008
 /  modified1: 
 /
-/  PURPOSE:Finds all shock mach numbers and injects energy into CR color 
-/  fields
+/  PURPOSE:Finds all shock mach numbers 
 /
 /  RETURNS:
 /    SUCCESS or FAIL
   if (NumberOfBaryonFields == 0)
     return SUCCESS;
  
-  this->DebugCheck("FindShocks");
-
-  long shockcounter, stopcounter;
+  this->DebugCheck((char *)"FindShocks");
 
   float invdx = 1./(CellWidth[0][0]);
   float inv2dx = 1./(2.*CellWidth[0][0]);
 
   int i, j, k, index,
     tempi, posti, prei;
-  float mu,
-    preT, postT, tempjumpmag,
+  float preT, postT, tempjumpmag,
     gradtx, gradty, gradtz,
     maxdiv, temprat, tempmach;
-  float sarea,vol,
-    Csound;
   float num;
-  double energyin;
-
-  float prepop,CRfraction;
-  int thisprepopbin,thismachbin;
 
   int DensNum, TENum, GENum, 
     Vel1Num, Vel2Num, Vel3Num;
 
-  int MachNum, CRNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
 
   /* Compute size (in floats) of the current grid. */
   int size = 1;
     ENZO_FAIL("Error in IdentifyPhysicalQuantities.");
   }
    
-  // Get CR species fields.
-
-  if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
-      ENZO_FAIL("Error in IdentifyCRSpeciesFields.");
+  // Get Shock species fields.
+  if (IdentifyShockSpeciesFields(MachNum,PSTempNum,PSDenNum) == FAIL) {
+      ENZO_FAIL("Error in IdentifyShockSpeciesFields.");
   }
 
   /* Get easy to handle pointers for each variable. */
  
   float *density     = BaryonField[DensNum];
-  float *totalenergy = BaryonField[TENum];
-  float *gasenergy   = BaryonField[GENum];
   float *velocity1   = BaryonField[Vel1Num];
   float *velocity2   = BaryonField[Vel2Num];
   float *velocity3   = BaryonField[Vel3Num];
   float *mach        = BaryonField[MachNum];
-  float *cr          = BaryonField[CRNum];
   float *pstemp      = BaryonField[PSTempNum];
   float *psden       = BaryonField[PSDenNum];
 
   }
   
   float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1,
-    VelocityUnits = 1, TimeUnits = 1, aUnits = 1;
+    VelocityUnits = 1, TimeUnits = 1;
 
   if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
 	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
   }
  
   // calculate cell temperature, set default values for mach
-  float energy;
   for (i=0;i<size;i++){
-    if(CRModel == 2)
-      cr[i] = 0.0;
-
     entropy[i] = tempgrad_dot_entropygrad[i] = mach[i] = 0.0;
     flowdivergence[i] = (double)(0.0);
     entropy[i] = temperature[i] / (pow(density[i],(Gamma - 1.0)));
     }
   }
 
-  //  fprintf(stdout, "ShockTemperatureFloor= %e TemperatureUnits= %e\n",ShockTemperatureFloor,
-  //  TemperatureUnits);
-
-
   //Calculate temperature gradient dotted with entropy gradient
   //Calculate the flow divergence.
   //Do this for all cells except last ghost zone.
       cell, and stop looking
    Make sure that the temperature/density keep increasing or decreasing
    After the ends are found, compute the mach number.
-   Then compute the CR energy injected, multiplied by the timestep!
 
    In order to get statistics in terms of pre-shock quantities later, put 
      the shock in the pre-shock cell(i.e. mach number, cr...)
 
 	mach[index] = tempmach;
 
-
-	//Calculate mass flux self-consistently
-	temprat = postT/preT;
-	tempmach = 
-	  sqrt(( 8.0*temprat - 7.0e0 + 
-		 sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat)
-		       + 15.0e0) )/5.0e0);
-	//Speed of sound in code units of velocity
-	Csound = sqrt(Gamma*kboltz*preT/(Mu*mh))/VelocityUnits;
-
-	/*-------------Lots of Extra Crap------------/
-	sarea = CellWidth[0][0] * CellWidth[0][0];
-	vol = CellWidth[0][0] * sarea;
-	//For now calculate the kinetic energy flux through the shock.
-
-	//This is in code units of energy.  
-	//Should multiply by vol*DensityUnits*VelUnits^2
-	//to get units of ergs
-	energyin = 0.5*1.19*sarea*density[prei]*
-	  pow(Csound*mach[index],3.0)*dtFixed;
-
-	//Now in honest to god ergs per cell
-	energyin *= LengthUnits*LengthUnits*DensityUnits*
-	  VelocityUnits*VelocityUnits*VelocityUnits*TimeUnits;
-
-	//Now convert to the same as gas_energy:
-	energyin /= vol*LengthUnits*LengthUnits*LengthUnits;
-	energyin /= density[index]*DensityUnits;
-	energyin /= VelocityUnits*VelocityUnits;
-	/--------------------------------------------*/
-	//*1.19e0
-	//Put energy in postshock cell
-	energyin = 0.5e0*pow(Csound*tempmach,3.0)*
-	  (density[prei]/density[posti])*dtFixed*
-	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-
-	//Put energy in middle cell
-// 	energyin = 0.5e0*pow(Csound*tempmach,3.0)*
-// 	  (density[prei]/density[index])*dtFixed*
-// 	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-	
-	prepop = cr[prei]/gasenergy[prei];
-	
-	thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)*
-			    (log10(tempmach) - CosmicRayData.CRMinMach)/
-			    (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach));
-	
-	thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)*
-			      (prepop - CosmicRayData.CRMinPrePop)/
-			      (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop));
-	
-	thismachbin = max(0, (int)(thismachbin));
-	thismachbin = min((CosmicRayData.CRNumberMachValues -1),
-			  (int)(thismachbin));			 
-	
-	thisprepopbin = max(0, (int)(thisprepopbin));
-	thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1),
-			  (int)(thisprepopbin));			 
-
-	CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin];
-
-	//Want this in Code Units
-	if(CRModel == 1)
-	  cr[posti]+= CRfraction*energyin;
-	//	  cr[index]+= CRfraction*energyin;
-
-	if(CRModel == 2)
-	  cr[index] = CRfraction*energyin/(dtFixed*TimeUnits);
-	
-	if(CRModel == 3){
-	  if(CRfraction*energyin > gasenergy[index] || 
-	     CRfraction*energyin > totalenergy[index]){
-	    gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    cr[index] += 0.9*min(gasenergy[index],totalenergy[index]);
-	  } else {
-	    cr[index] += CRfraction*energyin;
-	    gasenergy[index] -= CRfraction*energyin;
-	    totalenergy[index] -= CRfraction*energyin;
-	  }
-	}
 	if(StorePreShockFields){
 	  pstemp[index] = max(temperature[prei],ShockTemperatureFloor);
 	  psden[index] = density[prei];
     }
   }
   
-  /* deallocate temporary space for solver */
+  /* deallocate temporary space */
   
   delete [] temperature;
   delete [] flowdivergence;
   if (NumberOfBaryonFields == 0)
     return SUCCESS;
  
-  this->DebugCheck("FindShocks");
+  this->DebugCheck((char *)"FindShocks");
 
-  long shockcounter, stopcounter;
-
-  float invdx = 1./(CellWidth[0][0]);
   float inv2dx = 1./(2.*CellWidth[0][0]);
-  float inv2dx2 = invdx*invdx;
 
   int is=GridStartIndex[0];
   int js=GridStartIndex[1];
 
   int i, j, k, index,
     tempi, posti, prei;
-  float mu,centervelx,centervely,centervelz,
+  float Csound, centervelx,centervely,centervelz,
     preV, postV,veljumpmag,
     v1jump,v2jump,v3jump,
     gradvx, gradvy, gradvz,
     maxdiv, thisjump, oldjump, velmach;
-  float sarea,vol,
-    Csound;
   float num;
-  double energyin;
-
-  float prepop,CRfraction;
-  int thisprepopbin,thismachbin;
 
   int DensNum, TENum, GENum, 
     Vel1Num, Vel2Num, Vel3Num;
 
-  int MachNum, CRNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
 
   /* Compute size (in floats) of the current grid. */
   int size = 1;
    
   // Get CR species fields.
 
-  if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
+  if (IdentifyShockSpeciesFields(MachNum,PSTempNum,PSDenNum) == FAIL) {
     ENZO_FAIL("Error in IdentifyCRSpeciesFields.");
   }
 
   /* Get easy to handle pointers for each variable. */
  
   float *density     = BaryonField[DensNum];
-  float *totalenergy = BaryonField[TENum];
-  float *gasenergy   = BaryonField[GENum];
   float *velocity1   = BaryonField[Vel1Num];
   float *velocity2   = BaryonField[Vel2Num];
   float *velocity3   = BaryonField[Vel3Num];
   float *mach        = BaryonField[MachNum];
-  float *cr          = BaryonField[CRNum];
   float *pstemp      = BaryonField[PSTempNum];
   float *psden       = BaryonField[PSDenNum];
 
   }
   
   float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1,
-    VelocityUnits = 1, TimeUnits = 1,  aUnits = 1;
+    VelocityUnits = 1, TimeUnits = 1;
 
   if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
 	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
   }
  
   // calculate cell temperature, set default values for mach
-  float energy;
   for (i=0;i<size;i++){
-    if(CRModel == 2)
-      cr[i] = 0.0;
-
     entropy[i] = tempgrad_dot_entropygrad[i] = mach[i] = 0.0;
     flowdivergence[i] = (double)(0.0);
     entropy[i] = temperature[i] / (pow(density[i],(Gamma - 1.0)));
 	    break;
 	  }
 	  //Make sure density/temperature keeps increasing
-// 	  if(temperature[posti] < temperature[tempi] ||
-// 	     density[posti] < density[tempi]){
-// 	    posti = tempi;
-// 	    break;
-// 	  }
+	  if(temperature[posti] < temperature[tempi] ||
+	     density[posti] < density[tempi]){
+	    posti = tempi;
+	    break;
+	  }
 	  oldjump = thisjump;
 	  //Check for a shock in the current cell.  If not, set postV
 	  //and break out. Use last actual shocked cell as cell, not
 	    break;
 	  }
 	  //Make sure density/temperature keeps decreasing
-// 	  if(temperature[prei] > temperature[tempi] ||
-// 	     density[prei] > density[tempi]){
-// 	    posti = tempi;
-// 	    break;
-// 	  }
+	  if(temperature[prei] > temperature[tempi] ||
+	     density[prei] > density[tempi]){
+	    posti = tempi;
+	    break;
+	  }
 	  //Check for a shock in the current cell.  If not, set preV
 	  //and break out.
 	  if(flowdivergence[prei] >= 0.0){
 	if(num == -1)
 	  continue;
 
-// 	temprat = max(ShockTemperatureFloor,postV)/(max(ShockTemperatureFloor,preV));
-// 	//temprat = max(postV,ShockTemperatureFloor)/(max(preV,ShockTemperatureFloor));
+	//temprat = max(ShockTemperatureFloor,postV)/(max(ShockTemperatureFloor,preV));
+	//temprat = max(postV,ShockTemperatureFloor)/(max(preV,ShockTemperatureFloor));
 	
-//  	if(max(temperature[posti],ShockTemperatureFloor) <= max(temperature[prei],ShockTemperatureFloor))
-//  	  continue;
+ 	if(max(temperature[posti],ShockTemperatureFloor) <= max(temperature[prei],ShockTemperatureFloor))
+ 	  continue;
 
-//  	if(density[posti] <= density[prei])
-//  	  continue;
+ 	if(density[posti] <= density[prei])
+ 	  continue;
 
 	v1jump = gradvx*(velocity1[posti] - velocity1[prei]);
 	if(GridRank > 1)
 
 	mach[index] = velmach; 
 
-	energyin = 0.5e0*pow(thisjump,3.0)*
-	  (density[prei]/density[index])*dtFixed*
-	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-	
-	//Energy Injection
-	energyin = 0.5e0*pow(Csound*mach[index],3.0)*
-	  (density[prei]/density[index])*dtFixed*
-	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-
-	prepop = cr[prei]/gasenergy[prei];
-	
-	thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)*
-			    (log10(velmach) - CosmicRayData.CRMinMach)/
-			    (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach));
-	
-	thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)*
-			      (prepop - CosmicRayData.CRMinPrePop)/
-			      (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop));
-	
-	thismachbin = max(0, (int)(thismachbin));
-	thismachbin = min((CosmicRayData.CRNumberMachValues -1),
-			  (int)(thismachbin));			 
-	
-	thisprepopbin = max(0, (int)(thisprepopbin));
-	thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1),
-			  (int)(thisprepopbin));			 
-
-	CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin];
-
-//  	fprintf(stdout,"thisprepopbin: %i thismachbin: %i energyin: %e s
-//  	CRfraction:
-//  	%e\n",thisprepopbin,thismachbin,energyin,CRfraction);
-//  	fflush(stdout);
-	
-	//Want this in Code Units
-	if(CRModel == 1)
-	  cr[index]+= CRfraction*energyin;
-
-	if(CRModel == 2)
-	  cr[index] = CRfraction*energyin/(dtFixed*TimeUnits);
-	
-	if(CRModel == 3){
-	  if(CRfraction*energyin > gasenergy[index] || 
-	     CRfraction*energyin > totalenergy[index]){
-	    gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    cr[index] += 0.9*min(gasenergy[index],totalenergy[index]);
-	  } else {
-	    cr[index] += CRfraction*energyin;
-	    gasenergy[index] -= CRfraction*energyin;
-	    totalenergy[index] -= CRfraction*energyin;
-	  }
-	}
 	if(StorePreShockFields){
 	  pstemp[index] = max(temperature[prei],ShockTemperatureFloor);
 	  psden[index] = density[prei];
     }
   }
   
-  /* deallocate temporary space for solver */
+  /* deallocate temporary space */
   
   delete [] temperature;
   delete [] flowdivergence;
   if (NumberOfBaryonFields == 0)
     return SUCCESS;
  
-  this->DebugCheck("FindShocks");
+  this->DebugCheck((char *)"FindShocks");
 
-  long shockcounter, stopcounter;
-
-  float invdx = 1./(CellWidth[0][0]);
   float inv2dx = 1./(2.*CellWidth[0][0]);
-  float inv2dx2 = invdx*invdx;
 
   int is=GridStartIndex[0];
   int js=GridStartIndex[1];
   int ke=GridEndIndex[2];
 
   int i, j, k, index,
-    tempi, posti, prei;
-  float mu,centervelx,centervely,centervelz,
-    preV, postV,veljumpmag,
-    v1jump,v2jump,v3jump,
-    gradvx, gradvy, gradvz,
-    maxdiv, thisjump, oldjump, velmach;
-  float sarea,vol,
-    Csound;
-  float num;
-  double energyin;
-
-  float prepop,CRfraction;
-  int thisprepopbin,thismachbin;
+    prei;
+  float v1jump, v2jump, v3jump, Csound;
 
   int DensNum, TENum, GENum, 
     Vel1Num, Vel2Num, Vel3Num;
 
-  int MachNum, CRNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
 
   /* Compute size (in floats) of the current grid. */
   int size = 1;
    
   // Get CR species fields.
 
-  if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
-    ENZO_FAIL("Error in IdentifyCRSpeciesFields.");
+  if (IdentifyShockSpeciesFields(MachNum,PSTempNum,PSDenNum) == FAIL) {
+    ENZO_FAIL("Error in IdentifyShockSpeciesFields.");
   }
 
   /* Get easy to handle pointers for each variable. */
  
   float *density     = BaryonField[DensNum];
-  float *totalenergy = BaryonField[TENum];
-  float *gasenergy   = BaryonField[GENum];
   float *velocity1   = BaryonField[Vel1Num];
   float *velocity2   = BaryonField[Vel2Num];
   float *velocity3   = BaryonField[Vel3Num];
   float *mach        = BaryonField[MachNum];
-  float *cr          = BaryonField[CRNum];
   float *pstemp      = BaryonField[PSTempNum];
   float *psden       = BaryonField[PSDenNum];
 
   }
   
   float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1,
-    VelocityUnits = 1, TimeUnits = 1, aUnits = 1;
+    VelocityUnits = 1, TimeUnits = 1;
 
   if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
 	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
   }
  
   // calculate cell temperature, set default values for mach
-  float energy;
   for (i=0;i<size;i++){
-    if(CRModel == 2)
-      cr[i] = 0.0;
-
     entropy[i] = tempgrad_dot_entropygrad[i] = mach[i] = 0.0;
     flowdivergence[i] = (double)(0.0);
     entropy[i] = temperature[i] / (pow(density[i],(Gamma - 1.0)));
 	  continue;
 
 	//x-direction
-
-	v1jump = fabs(velocity1[index]-velocity1[index-1]);
-	if(velocity1[index] > velocity1[index-1])
+	v1jump = fabs(velocity1[index+1]-velocity1[index-1]);
+	if(velocity1[index+1] > velocity1[index-1])
 	  prei = index-1;
 	else
-	  prei = index;
+	  prei = index+1;
 	
 	Csound = sqrt(Gamma*kboltz*temperature[prei]/
 		      (Mu*mh))/VelocityUnits;
 	mach[prei] = (1.0/3.0)*(2.0*v1jump/Csound + 
 				sqrt(9.0+4.0*v1jump*v1jump/Csound/Csound));
-      	if(flowdivergence[prei] >= 0.0)
-	  mach[prei]=0.0;
 	
 	//y-direction
 
 	if (GridRank > 1){
-	  v2jump = fabs(velocity2[index]-velocity2[index-GridDimension[0]]);
-	  if(velocity2[index] > velocity2[index-GridDimension[0]])
+	  v2jump = fabs(velocity2[index+GridDimension[0]]-velocity2[index-GridDimension[0]]);
+	  if(velocity2[index+GridDimension[0]] > velocity2[index-GridDimension[0]])
 	    prei = index-GridDimension[0];
 	  else
-	    prei = index;
+	    prei = index+GridDimension[0];
 	  Csound = sqrt(Gamma*kboltz*temperature[prei]/
 			(Mu*mh))/VelocityUnits;
 	  mach[prei] = sqrt(mach[prei]*mach[prei] + 
 				       sqrt(9.0+4.0*v2jump*v2jump/Csound/Csound))*
 			    (1.0/3.0)*(2.0*v2jump/Csound + 
 				       sqrt(9.0+4.0*v2jump*v2jump/Csound/Csound)));
-	  if(flowdivergence[prei] >= 0.0)
-	    mach[prei]=0.0;
 	}	
 	
 	//z-direction
 
 	if (GridRank > 2){
-	  v3jump = fabs(velocity3[index]-velocity3[index-GridDimension[0]]);
-	  if(velocity3[index] > velocity3[index-
-					  GridDimension[0]*GridDimension[1]])
+	  v3jump = fabs(velocity3[index+GridDimension[0]*GridDimension[1]]-
+			velocity3[index-GridDimension[0]*GridDimension[1]]);
+	  if(velocity3[index+GridDimension[0]*GridDimension[1]] > 
+	     velocity3[index-GridDimension[0]*GridDimension[1]])
 	    prei = index-GridDimension[0]*GridDimension[1];
 	  else
-	    prei = index;
+	    prei = index+GridDimension[0]*GridDimension[1];
 	  Csound = sqrt(Gamma*kboltz*temperature[prei]/
 			(Mu*mh))/VelocityUnits;
 	  mach[prei] = sqrt(mach[prei]*mach[prei] + 
 				       sqrt(9.0+4.0*v3jump*v3jump/Csound/Csound))*
 			    (1.0/3.0)*(2.0*v3jump/Csound + 
 				       sqrt(9.0+4.0*v3jump*v3jump/Csound/Csound)));
-	  if(flowdivergence[prei] >= 0.0)
-	    mach[prei]=0.0;
 	}	
-      }
-    }
-  }
-  for(k=ks; k<=ke;k++){
-    for(j=js; j<=je;j++){
-      for(i=is; i<=ie;i++){
-	
- 	index = i + GridDimension[0]*(j + GridDimension[1]*k);	
-	Csound = sqrt(Gamma*kboltz*temperature[index]/
-		      (Mu*mh))/VelocityUnits;	
-	//Energy Injection
-	energyin = 0.5e0*pow(Csound*mach[index],3.0)*
-	  (density[index]/density[index])*dtFixed*
-	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-	
-	prepop = cr[index]/gasenergy[index];
-	
-	thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)*
-			    (log10(mach[index]) - CosmicRayData.CRMinMach)/
-			    (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach));
-	
-	thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)*
-			      (prepop - CosmicRayData.CRMinPrePop)/
-			      (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop));
-	
-	thismachbin = max(0, (int)(thismachbin));
-	thismachbin = min((CosmicRayData.CRNumberMachValues -1),
-			  (int)(thismachbin));			 
-	
-	thisprepopbin = max(0, (int)(thisprepopbin));
-	thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1),
-			    (int)(thisprepopbin));			 
-	
-	CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin];
-	
-	//Want this in Code Units
-	if(CRModel == 1)
-	  cr[index]+= CRfraction*energyin;	
-	
-	if(CRModel == 2)
-	  cr[index] = CRfraction*energyin/(dtFixed*TimeUnits);
-	
-	if(CRModel == 3){
-	  if(CRfraction*energyin > gasenergy[index] || 
-	     CRfraction*energyin > totalenergy[index]){
-	    gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    cr[index] += 0.9*min(gasenergy[index],totalenergy[index]);
-	  } else {
-	    cr[index] += CRfraction*energyin;
-	    gasenergy[index] -= CRfraction*energyin;
-	    totalenergy[index] -= CRfraction*energyin;
-	  }
-	}
+
 	if(StorePreShockFields){
 	  pstemp[index] = max(temperature[prei],ShockTemperatureFloor);
 	  psden[index] = density[prei];
 	}	
+	
       }
     }
   }
   
-  /* deallocate temporary space for solver */
+  /* deallocate temporary space */
   
   delete [] temperature;
   delete [] flowdivergence;
   if (NumberOfBaryonFields == 0)
     return SUCCESS;
  
-  this->DebugCheck("FindShocks");
-
-  long shockcounter, stopcounter;
+  this->DebugCheck((char *)"FindShocks");
 
   float invdx = 1./(CellWidth[0][0]);
   float inv2dx = 1./(2.*CellWidth[0][0]);
   int je=GridEndIndex[1];
   int ke=GridEndIndex[2];
 
-  int i, j, k, index,
-    tempi, posti, prei;
-  float mu,centervelx,centervely,centervelz,
-    preV, postV,veljumpmag,
-    v1jump,v2jump,v3jump,
-    gradvx, gradvy, gradvz,
-    maxdiv, thisjump, oldjump, velmach;
-  float sarea,vol,
-    Csound;
-  float num;
-  double energyin;
-
-  float prepop,CRfraction;
-  int thisprepopbin,thismachbin;
+  int i, j, k, index,posti, prei;
+  float maxdiv;
 
   //Specific for Split Temperature:
   float tempden, temptemp, mach1,
-    postden, preden, postT, preT,
+    postT, preT,
     temprat;
   int centerfound;
 
   int DensNum, TENum, GENum, 
     Vel1Num, Vel2Num, Vel3Num;
 
-  int MachNum, CRNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
 
   /* Compute size (in floats) of the current grid. */
   int size = 1;
    
   // Get CR species fields.
 
-  if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
-    ENZO_FAIL("Error in IdentifyCRSpeciesFields.");
+  if (IdentifyShockSpeciesFields(MachNum,PSTempNum,PSDenNum) == FAIL) {
+    ENZO_FAIL("Error in IdentifyShockSpeciesFields.");
   }
 
   /* Get easy to handle pointers for each variable. */
  
   float *density     = BaryonField[DensNum];
-  float *totalenergy = BaryonField[TENum];
-  float *gasenergy   = BaryonField[GENum];
   float *velocity1   = BaryonField[Vel1Num];
   float *velocity2   = BaryonField[Vel2Num];
   float *velocity3   = BaryonField[Vel3Num];
   float *mach        = BaryonField[MachNum];
-  float *cr          = BaryonField[CRNum];
   float *pstemp      = BaryonField[PSTempNum];
   float *psden       = BaryonField[PSDenNum];
 
   }
   
   float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1,
-    VelocityUnits = 1, TimeUnits = 1, aUnits = 1;
+    VelocityUnits = 1, TimeUnits = 1;
 
   if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
 	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
   }
  
   // calculate cell temperature, set default values for mach
-  float energy;
   for (i=0;i<size;i++){
-    if(CRModel == 2)
-      cr[i] = 0.0;
-
     entropy[i] = tempgrad_dot_entropygrad[i] = mach[i] = 0.0;
     flowdivergence[i] = (double)(0.0);
     entropy[i] = temperature[i] / (pow(density[i],(Gamma - 1.0)));
 
 	//x-direction
 
+	posti = index;
+	tempden = density[index];
+	maxdiv = flowdivergence[index];
+	temptemp = temperature[index];
 	while(true){
-	  posti = index;
-	  tempden = density[index];
-	  maxdiv = flowdivergence[index];
-	  temptemp = temperature[index];
 	  if(temperature[index+1] >= temperature[index-1]){
 	    posti++;
 	  }else{
 	  }
 	  maxdiv = flowdivergence[posti];
 	}
+	prei = index;
+	tempden = density[index];
+	temptemp = temperature[index];
 	while(true){
-	  prei = index;
-	  tempden = density[index];
-	  temptemp = temperature[index];
 	  if(temperature[index+1] < temperature[index-1]){
 	    prei++;
 	  }else{
 	
 	//y-direction
 	if(GridRank > 1){
+	  posti = index;
+	  tempden = density[index];
+	  temptemp = temperature[index];
 	  while(true){
-	    posti = index;
-	    tempden = density[index];
-	    temptemp = temperature[index];
 	    maxdiv = flowdivergence[index];
 	    if(temperature[index+GridDimension[0]] >= 
 	       temperature[index-GridDimension[0]]){
 	    }
 	    maxdiv = flowdivergence[posti];
 	  }
+	  prei = index;
+	  tempden = density[index];
+	  temptemp = temperature[index];
 	  while(true){
-	    prei = index;
-	    tempden = density[index];
-	    temptemp = temperature[index];
 	    if(temperature[index+GridDimension[0]] < 
 	       temperature[index-GridDimension[0]]){
 	      prei+=GridDimension[0];
 	
 	//z-direction
 	if(GridRank > 2){
+	  posti = index;
+	  tempden = density[index];
+	  temptemp = temperature[index];
 	  while(true){
-	    posti = index;
-	    tempden = density[index];
-	    temptemp = temperature[index];
 	    maxdiv = flowdivergence[index];
 	    if(temperature[index+GridDimension[0]*GridDimension[1]] >= 
 	       temperature[index-GridDimension[0]*GridDimension[1]]){
 	    }
 	    maxdiv = flowdivergence[posti];
 	  }
+	  prei = index;
+	  tempden = density[index];
+	  temptemp = temperature[index];
 	  while(true){
-	    prei = index;
-	    tempden = density[index];
-	    temptemp = temperature[index];
 	    if(temperature[index+GridDimension[0]*GridDimension[1]] < 
 	       temperature[index-GridDimension[0]*GridDimension[1]]){
 	      prei+=GridDimension[0]*GridDimension[1];
 	      mach[index]=mach1;
 	  }
 	}  // GridRank > 2
-      } // for i
-    } // for j
-  } // for k
-  for(k=ks; k<=ke;k++){
-    for(j=js; j<=je;j++){
-      for(i=is; i<=ie;i++){
-	
- 	index = i + GridDimension[0]*(j + GridDimension[1]*k);	
-	Csound = sqrt(Gamma*kboltz*temperature[index]/
-		      (Mu*mh))/VelocityUnits;	
-	//Energy Injection
-	energyin = 0.5e0*pow(Csound*mach[index],3.0)*
-	  (density[index]/density[index])*dtFixed*
-	  TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0];
-	
-	prepop = cr[index]/gasenergy[index];
-	
-	thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)*
-			    (log10(mach[index]) - CosmicRayData.CRMinMach)/
-			    (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach));
-	
-	thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)*
-			      (prepop - CosmicRayData.CRMinPrePop)/
-			      (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop));
-	
-	thismachbin = max(0, (int)(thismachbin));
-	thismachbin = min((CosmicRayData.CRNumberMachValues -1),
-			  (int)(thismachbin));			 
-	
-	thisprepopbin = max(0, (int)(thisprepopbin));
-	thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1),
-			    (int)(thisprepopbin));			 
-	
-	CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin];
-	
-	//Want this in Code Units
-	if(CRModel == 1)
-	  cr[index]+= CRfraction*energyin;	
-	
-	if(CRModel == 2)
-	  cr[index] = CRfraction*energyin/(dtFixed*TimeUnits);
-	
-	if(CRModel == 3){
-	  if(CRfraction*energyin > gasenergy[index] || 
-	     CRfraction*energyin > totalenergy[index]){
-	    gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]);
-	    cr[index] += 0.9*min(gasenergy[index],totalenergy[index]);
-	  } else {
-	    cr[index] += CRfraction*energyin;
-	    gasenergy[index] -= CRfraction*energyin;
-	    totalenergy[index] -= CRfraction*energyin;
-	  }
-	}
 	if(StorePreShockFields){
 	  pstemp[index] = max(temperature[prei],ShockTemperatureFloor);
 	  psden[index] = density[prei];
 	}
-      }
-    }
-  }
+      } // for i
+    } // for j
+  } // for k
   
-  /* deallocate temporary space for solver */
+  /* deallocate temporary space */
   
   delete [] temperature;
   delete [] flowdivergence;
   delete [] tempgrad_dot_entropygrad;
   delete [] entropy;
-  
+
   return SUCCESS;
 }

File src/enzo/Grid_IdentifyCRSpeciesFields.C

-/***********************************************************************
-/
-/  GRID CLASS (IDENTIFY THE SPECIES FIELDS FOR SAM SKILLMAN'S COSMIC RAYS)
-/
-/  written by: Samuel Skillman
-/  date:       May, 2008
-/  modified1:
-/
-/  PURPOSE:
-/
-/  NOTE:
-/
-************************************************************************/
- 
-#include <stdio.h>
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "fortran.def"
- 
-/* function prototypes */
- 
-int FindField(int f, int farray[], int n);
- 
-int grid::IdentifyCRSpeciesFields(int &MachNum,int &CRNum, 
-				  int &PSTempNum, int &PSDenNum)
-{
-  MachNum = CRNum = PSTempNum = PSDenNum = 0;
-
-  // Basic: Mach, CR Protons
-  if (CRModel) {
-    MachNum  = FindField(Mach, FieldType, NumberOfBaryonFields);
-    CRNum    = FindField(CRDensity, FieldType, NumberOfBaryonFields);
-    if (StorePreShockFields){
-      PSTempNum= FindField(PreShockTemperature, FieldType, NumberOfBaryonFields);
-      PSDenNum = FindField(PreShockDensity, FieldType, NumberOfBaryonFields);    
-    }     
-  }
-
-  if ((MachNum < 0) || (CRNum < 0) || (PSTempNum < 0) || (PSDenNum < 0)) {
-    fprintf(stderr,"Error identifying species for CRModel = %"ISYM" MachNum= %"ISYM" CRNum = %"ISYM" PSTempNum = %"ISYM" PSDenNum = %"ISYM" NBaryonFs = %"ISYM".\n",
-	    CRModel,MachNum,CRNum,PSTempNum,PSDenNum,NumberOfBaryonFields);
-    return FAIL;
-  }
-  return SUCCESS;
-}

File src/enzo/Grid_IdentifyShockSpeciesFields.C

+/***********************************************************************
+/
+/  GRID CLASS (IDENTIFY THE SHOCK SPECIES FIELDS)
+/
+/  written by: Samuel Skillman
+/  date:       May, 2008
+/  modified1:
+/
+/  PURPOSE:
+/
+/  NOTE:
+/
+************************************************************************/
+ 
+#include <stdio.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "fortran.def"
+ 
+/* function prototypes */
+ 
+int FindField(int f, int farray[], int n);
+ 
+int grid::IdentifyShockSpeciesFields(int &MachNum, int &PSTempNum, int &PSDenNum)
+{
+  MachNum = PSTempNum = PSDenNum = 0;
+
+  // Basic: Mach, CR Protons
+  if (ShockMethod) {
+    MachNum  = FindField(Mach, FieldType, NumberOfBaryonFields);
+    if (StorePreShockFields){
+      PSTempNum= FindField(PreShockTemperature, FieldType, NumberOfBaryonFields);
+      PSDenNum = FindField(PreShockDensity, FieldType, NumberOfBaryonFields);    
+    }     
+  }
+
+  if ((MachNum < 0) || (PSTempNum < 0) || (PSDenNum < 0)) {
+    fprintf(stderr,"Error identifying species for ShockMethod = %"ISYM" MachNum= %"ISYM" PSTempNum = %"ISYM" PSDenNum = %"ISYM" NBaryonFs = %"ISYM".\n",
+	    ShockMethod,MachNum,PSTempNum,PSDenNum,NumberOfBaryonFields);
+    return FAIL;
+  }
+  return SUCCESS;
+}

File src/enzo/Grid_InterpolateBoundaryFromParent.C

  
       if (FieldType[field] == Density)
 	FieldPointer = TemporaryDensityField;
-      else if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
+      else 
 	FieldPointer = TemporaryField;
  
       /* Copy needed portion of temp field to current grid. */

File src/enzo/Grid_InterpolateFieldValues.C

   float *TemporaryField, *TemporaryDensityField, *Work,
         *ParentTemp[MAX_NUMBER_OF_BARYON_FIELDS], *FieldPointer;
  
-  int MyInterpolationMethod = InterpolationMethod;   
   if (NumberOfBaryonFields > 0) {
 
     interp_error = FALSE;
        is done for the entire current grid, not just it's boundaries.
        (skip density since we did it already) */
  
-      if (FieldTypeNoInterpolate(FieldType[field]) == FALSE){
-        if (HydroMethod == Zeus_Hydro){
-	        InterpolationMethod = (SecondOrderBFlag[field] == 0) ?
-	            SecondOrderA : SecondOrderC;
-        }
-      } else {
-        /* Use nearest grid point interpolation for fields that 
-           shouldn't ever be averaged. */ 
-        MyInterpolationMethod = FirstOrderA; 
+      if (HydroMethod == Zeus_Hydro){
+	InterpolationMethod = (SecondOrderBFlag[field] == 0) ?
+	  SecondOrderA : SecondOrderC;
       }
       //      fprintf(stdout, "grid:: InterpolateBoundaryFromParent[4], field = %d\n", field); 
 

File src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C

   int iTE = ietot;
   int ExtraField[2];
   int ForbidNum;
-  int CRNum, MachNum, PSTempNum, PSDenNum;
+  int MachNum, PSTempNum, PSDenNum;
  
   inits_type *tempbuffer = NULL;
   int *int_tempbuffer = NULL;
       fprintf(stderr, "Initializing Forbidden Refinement color field\n");
       FieldType[ForbidNum = NumberOfBaryonFields++] = ForbiddenRefinement;
     }
-    if(CRModel){
+    if(ShockMethod){
       FieldType[MachNum   = NumberOfBaryonFields++] = Mach;
       if(StorePreShockFields){
 	FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature;
 	FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity;
       }
-      FieldType[CRNum     = NumberOfBaryonFields++] = CRDensity;
     }    
   }
 
 	      }
 	    }
 	}
+      } // end: if (CosmologySimulationDensityName != NULL)
 
-	  
-
-
-	//Shock/Cosmic Ray Model
-	if (CRModel && ReadData) {
+      // Shock/Cosmic Ray Model
+      if (ShockMethod && ReadData) {
+	for (i = 0; i < size; i++) {
 	  BaryonField[MachNum][i] = tiny_number;
-	  BaryonField[CRNum][i] = tiny_number;
 	  if (StorePreShockFields) {
 	    BaryonField[PSTempNum][i] = tiny_number;
 	    BaryonField[PSDenNum][i] = tiny_number;
-	  }
-	} // end: if (CRModel && ReadData)
- 
-      } // end: if (CosmologySimulationDensityName != NULL)
+	    }
+	}
+      } // end: if (ShockMethod && ReadData)
 
     } // end: if (NumberOfBaryonFields > 0)
 

File src/enzo/Grid_PhotonSortLinkedLists.C

+/***********************************************************************
+/
+/  GRID CLASS (SORT LINKED LISTS OF PHOTONS)
+/
+/  written by: Stephen Skory
+/  date:       June, 2011
+/  modified1:
+/
+/  PURPOSE:  Sorts the linked lists of arrays.
+/
+************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+
+Eint32 compare_ss(const void *a, const void *b);
+PhotonPackageEntry *LinkedListToArray(PhotonPackageEntry *Node, int n);
+PhotonPackageEntry* DeletePhotonPackage(PhotonPackageEntry *PP);
+void InsertPhotonAfter(PhotonPackageEntry * &Node, PhotonPackageEntry * &NewNode);
+
+int grid::PhotonSortLinkedLists(void)
+{
+  // Photons of all three types are sorted following the method used
+  // in Grid_MergePausedPhotonPackages.
+
+  if (MyProcessorNumber != ProcessorNumber) return SUCCESS;
+  
+  int nphotons, dim, count;
+  PhotonPackageEntry *PP, *TempPP, *NewPack;
+  
+  // PhotonPackages
+  nphotons = 0;
+  PP = PhotonPackages->NextPackage;
+  while (PP != NULL) {
+    PP = PP->NextPackage;
+    nphotons++;
+  }
+  
+  PP = PhotonPackages->NextPackage;
+  TempPP = LinkedListToArray(PP, nphotons);
+  // Sort the list.
+  qsort(TempPP, nphotons, sizeof(PhotonPackageEntry), compare_ss);
+  
+  // Now we need to clear out PhotonPackages before we fill it back up again.
+  for (PP = PhotonPackages->NextPackage; PP; PP = PP->NextPackage)
+    PP = DeletePhotonPackage(PP);
+  PhotonPackages->NextPackage = NULL;
+  PhotonPackages->PreviousPackage = NULL;
+
+  // Fill it back in from the sorted list.
+  PP = PhotonPackages;
+  for (count = 0; count < nphotons; count++) {
+    NewPack = new PhotonPackageEntry;
+    NewPack->CurrentSource = TempPP[count].CurrentSource;
+    NewPack->Photons = TempPP[count].Photons;
+    NewPack->Type = TempPP[count].Type;
+    NewPack->Energy = TempPP[count].Energy;
+    NewPack->CrossSection = TempPP[count].CrossSection; //
+    NewPack->EmissionTimeInterval = TempPP[count].EmissionTimeInterval;
+    NewPack->EmissionTime = TempPP[count].EmissionTime;
+    NewPack->CurrentTime = TempPP[count].CurrentTime;
+    NewPack->Radius = TempPP[count].Radius;
+    NewPack->ColumnDensity = TempPP[count].ColumnDensity;
+    NewPack->ipix = TempPP[count].ipix;
+    NewPack->level = TempPP[count].level;
+    NewPack->SourcePositionDiff = TempPP[count].SourcePositionDiff;
+    for (dim = 0; dim < 3; dim++) {
+      NewPack->SourcePosition[dim] = TempPP[count].SourcePosition[dim];
+    }
+    InsertPhotonAfter(PP, NewPack);
+    PP = PP->NextPackage;
+  }
+  
+  // clean up
+  delete [] TempPP;
+  
+
+  // FinishedPhotonPackages
+  nphotons = 0;
+  PP = FinishedPhotonPackages->NextPackage;
+  while (PP != NULL) {
+    PP = PP->NextPackage;
+    nphotons++;
+  }
+  
+  PP = FinishedPhotonPackages->NextPackage;
+  TempPP = LinkedListToArray(PP, nphotons);
+  // Sort the list.
+  qsort(TempPP, nphotons, sizeof(PhotonPackageEntry), compare_ss);
+  
+  // Now we need to clear out FinishedPhotonPackages
+  // before we fill it back up again.
+  for (PP = FinishedPhotonPackages->NextPackage; PP; PP = PP->NextPackage)
+    PP = DeletePhotonPackage(PP);
+  FinishedPhotonPackages->NextPackage = NULL;
+  FinishedPhotonPackages->PreviousPackage = NULL;
+
+  // Fill it back in from the sorted list.
+  PP = FinishedPhotonPackages;
+  for (count = 0; count < nphotons; count++) {