Commits

John Wise  committed 599f469 Merge

Merging

  • Participants
  • Parent commits 51557cf, 0578407
  • Branches week-of-code

Comments (0)

Files changed (13)

File src/enzo/EvolvePhotons.C

 	  Temp->GridData->FinalizeRadiationFields();
     END_PERF(8);
 
-    START_PERF();
-    for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
-      for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel)
-	if (Temp->GridData->RadiationPresent() == TRUE) {
-
-	  if (RadiativeTransferCoupledRateSolver && 
-	      RadiativeTransferOpticallyThinH2) {
-	    Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
-	  } // ENDIF coupled & thinH2
-
-	  if (RadiativeTransferCoupledRateSolver) {
-	    int RTCoupledSolverIntermediateStep = TRUE;
-	    Temp->GridData->SolveRateAndCoolEquations(RTCoupledSolverIntermediateStep);
-	  }
-
-	  if (RadiativeTransferCoupledRateSolver &&
-	      RadiativeTransferInterpolateField)
-	    Temp->GridData->DeleteInterpolatedFields();
-
-	} /* ENDIF radiation */
-    END_PERF(9);
-
-    /* For the non-coupled (i.e. cells without radiation) rate & energy
-       solver, we have to set the H2 dissociation rates */
+    /* Set the optically-thin H2 dissociation rates */
 
     START_PERF();
     if (RadiativeTransferOpticallyThinH2)
       for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
 	for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel)
-	  if (Temp->GridData->RadiationPresent() == FALSE)
-	    Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
+	  Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
+    if (RadiativeTransferOpticallyThinXray)
+      for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
+	for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel)
+	  Temp->GridData->AddOpticallyThinXrays(AllStars, NumberOfSources);
     END_PERF(10);
 
+    START_PERF();
+    int RTCoupledSolverIntermediateStep = TRUE;
+    if (RadiativeTransferCoupledRateSolver)
+      for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
+	for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel)
+	  if (Temp->GridData->RadiationPresent() == TRUE) {
+
+	    Temp->GridData->SolveRateAndCoolEquations(RTCoupledSolverIntermediateStep);
+
+	  } /* ENDIF radiation */
+    END_PERF(9);
+
     /* Clean up temperature field */
 
     if (RadiationXRayComptonHeating)

File src/enzo/Grid_AddOpticallyThinXrays.C

+/***********************************************************************
+/
+/  (WRAPPER) ADD X-RAY EMISSION FROM RADIATIVE PARTICLES
+/
+/  written by: John Wise
+/  date:       March, 2012
+/  modified1:
+/
+/  PURPOSE:
+/
+************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.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 "CosmologyParameters.h"
+#include "Star.h"
+
+int grid::AddOpticallyThinXrays(Star *AllStars, int NumberOfSources)
+{
+
+  /* If we're merging rays, we already have a binary tree of the
+     sources.  We can use that to speed up the calculation when we
+     have more than 10 sources.  With smaller numbers, the overhead
+     makes the direct calculation faster. */
+
+#ifdef UNUSED // Difficult with multiple X-ray energy bins
+  if (RadiativeTransferSourceClustering == TRUE && NumberOfSources >= 10)
+    this->AddXraysFromTree();
+  else
+#endif
+  this->AddXraysFromSources(AllStars);
+
+  HasRadiation = TRUE;
+
+  return SUCCESS;
+
+}

File src/enzo/Grid_AddXraysFromSources.C

+/***********************************************************************
+/
+/  ADD X-RAYS EMISSION FROM SHINING PARTICLES
+/
+/  written by: John Wise
+/  date:       April, 2012
+/  modified1:
+/
+/  PURPOSE:
+/
+************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.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 "CosmologyParameters.h"
+#include "Star.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+FLOAT FindCrossSection(int type, float energy);
+
+int grid::AddXraysFromSources(Star *AllStars)
+{
+
+  const float EnergyThresholds[] = {13.6, 24.6, 54.4};
+  const float PopulationFractions[] = {1.0, 0.25, 0.25};
+  const float erg_eV = 1.602176e-12;
+  const double clight = 2.99792e10;
+  const double pc = 3.086e18;
+
+  Star *cstar;
+  FLOAT DomainWidth[MAX_DIMENSION];
+  FLOAT *ddr2[MAX_DIMENSION];
+  double Luminosity[MAX_ENERGY_BINS];
+  float energies[MAX_ENERGY_BINS], kph_r2, gamma_r2, dkph, dE;
+  int ipart, dim, a, i, j, k, bin, index, indixe, nbins;
+  int ActiveDims[MAX_DIMENSION];
+  int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
+      DINum, DIINum, HDINum;
+  double XrayLuminosity, LConv, CrossSectionConv;
+
+  if (MyProcessorNumber != ProcessorNumber)
+    return SUCCESS;
+
+  /* Exit if no star particles and not Photon Test */
+
+  if (AllStars == NULL && ProblemType != 50)
+    return SUCCESS;
+
+  /* Find Multi-species fields. */
+
+  if (this->IdentifySpeciesFields(DeNum, HINum, HIINum, HeINum, HeIINum, 
+				  HeIIINum, HMNum, H2INum, H2IINum, DINum, 
+				  DIINum, HDINum) == FAIL) {
+    ENZO_FAIL("Error in grid->IdentifySpeciesFields.\n");
+  }
+
+  /* Get photo-ionization fields */
+
+  int kphHINum, kphHeINum, kphHeIINum, kdissH2INum;
+  int gammaNum;
+  IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, kphHeIINum, 
+				  kdissH2INum);
+  const int kphNum[] = {kphHINum, kphHeINum, kphHeIINum};
+
+  /* If using cosmology, get units. */
+
+  float TemperatureUnits, DensityUnits, LengthUnits, VelocityUnits, 
+    TimeUnits, aUnits = 1;
+
+  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+	   &TimeUnits, &VelocityUnits, PhotonTime);
+
+  // Absorb the unit conversions into the cross-section
+  CrossSectionConv = (double)TimeUnits / ((double)LengthUnits * (double)LengthUnits);
+
+  // Convert from #/s to RT units
+  LConv = (double) TimeUnits / pow(LengthUnits,3);
+
+  for (dim = 0; dim < GridRank; dim++) {
+    DomainWidth[dim] = DomainRightEdge[dim] - DomainLeftEdge[dim];
+    ActiveDims[dim] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
+    ddr2[dim] = new FLOAT[ActiveDims[dim]];
+  }
+
+  /* Loop over radiation sources or star particles in the grid */
+
+  FLOAT CrossSections[3]; // HI, HeI, HeII
+  float xx, heat_factor, ion2_factor[3];
+  float nSecondaryHII = 1, nSecondaryHeII = 1;
+  double r2_inv;
+
+  ion2_factor[2] = 1.0;
+  if (RadiationXRaySecondaryIon == FALSE) {
+    heat_factor = 1.0;
+    ion2_factor[0] = ion2_factor[1] = 1.0;
+  }
+
+  if (ProblemType == 50) {
+
+    RadiationSourceEntry *RS;
+    for (RS = GlobalRadiationSources->NextSource; RS; RS = RS->NextSource) {
+
+      if (PhotonTime < RS->CreationTime && 
+	  PhotonTime > RS->CreationTime + RS->LifeTime)
+	continue;
+
+      /* Loop over energy bins and consider X-rays to be E >= 100 eV */
+
+      for (bin = 0; bin < RS->EnergyBins; bin++) {
+
+	if (RS->Energy[bin] < 100) continue;
+
+	XrayLuminosity = RS->SED[bin] * RS->Luminosity / LConv;
+	for (i = 0; i < 3; i++)
+	  CrossSections[i] = FindCrossSection(i, RS->Energy[bin]) * CrossSectionConv;
+
+	if (RadiationXRaySecondaryIon == TRUE) {
+	  nSecondaryHII = RS->Energy[bin] / 13.6;
+	  nSecondaryHeII = RS->Energy[bin] / 24.6;
+	}
+
+	/* Pre-calculate distances from cells to source */
+
+	for (dim = 0; dim < GridRank; dim++)
+	  for (i = 0, index = GridStartIndex[dim]; i < ActiveDims[dim]; 
+	       i++, index++) {
+
+	    // Calculate dr_i first, then square it
+	    ddr2[dim][i] = 
+	      fabs(CellLeftEdge[dim][index] + 0.5*CellWidth[dim][index] - 
+		   RS->Position[dim]);
+	    ddr2[dim][i] = min(ddr2[dim][i], DomainWidth[dim]-ddr2[dim][i]);
+	    ddr2[dim][i] = ddr2[dim][i] * ddr2[dim][i];
+	  }
+
+	/* Loop over absorbers then cells */
+
+	double radius2_yz;
+
+	for (k = 0; k < ActiveDims[2]; k++) {
+	  for (j = 0; j < ActiveDims[1]; j++) {
+	    radius2_yz = ddr2[1][j] + ddr2[2][k];
+	    index = GRIDINDEX(0, j, k);
+	    for (i = 0; i < ActiveDims[0]; i++, index++) {
+	      r2_inv = 1.0 / (radius2_yz + ddr2[0][i]);
+	      if (RadiationXRaySecondaryIon == TRUE) {
+		xx = max(BaryonField[HIINum][index] / 
+			 (BaryonField[HINum][index] + BaryonField[HIINum][index]),
+			 1e-4);
+		heat_factor = 0.9971 * (1.0 - powf(1.0 - powf(xx, 0.2663f), 1.3163));
+		ion2_factor[0] = 0.3908 * nSecondaryHII * 
+		  powf(1 - powf(xx, 0.4092f), 1.7592f);
+		ion2_factor[1] = 0.0554 * nSecondaryHeII * 
+		  powf(1 - powf(xx, 0.4614f), 1.6660f);
+	      }
+	      for (a = 0; a < 3; a++) {
+		dkph = (float) (XrayLuminosity * CrossSections[a] / (4.0 * M_PI)) * r2_inv;
+		dE = (RS->Energy[bin] - EnergyThresholds[a]);
+		BaryonField[kphNum[a]][index] += dkph * ion2_factor[a];
+		BaryonField[gammaNum][index] += dkph * dE * heat_factor;
+	      } // END: absorbers (a)
+	    } // END: i-direction
+	  } // END: j-direction
+	} // END: k-direction
+
+      } // ENDFOR bin
+      
+    } // ENDFOR sources
+
+  } // ENDIF ProblemType == 50
+
+  else {
+
+    for (cstar = AllStars; cstar; cstar = cstar->NextStar) {
+
+      // Skip if not 'living'
+      if (!(cstar->FeedbackFlag == NO_FEEDBACK ||
+	    cstar->FeedbackFlag == CONT_SUPERNOVA)) 
+	continue;
+      
+      /* Get energy bins and SED */
+
+      if (cstar->ComputePhotonRates(nbins, energies, Luminosity) == FAIL) {
+	ENZO_FAIL("Error in ComputePhotonRates.\n");
+      }
+
+      /* Pre-calculate distances from cells to source */
+
+      for (dim = 0; dim < GridRank; dim++)
+	for (i = 0, index = GridStartIndex[dim]; i < ActiveDims[dim]; 
+	     i++, index++) {
+	  
+	  // Calculate dr_i first, then square it
+	  ddr2[dim][i] = 
+	    fabs(CellLeftEdge[dim][index] + 0.5*CellWidth[dim][index] -
+		 cstar->pos[dim]);
+	  ddr2[dim][i] = min(ddr2[dim][i], DomainWidth[dim]-ddr2[dim][i]);
+	  ddr2[dim][i] = ddr2[dim][i] * ddr2[dim][i];
+	}
+
+      /* Loop over energy bins and consider X-rays to be E >= 100 eV */
+
+      for (bin = 0; bin < nbins; bin++) {
+
+	if (energies[bin] < 100) continue;
+
+	XrayLuminosity = Luminosity[bin] / LConv;
+	for (i = 0; i < 3; i++)
+	  CrossSections[i] = FindCrossSection(i, energies[bin]) * CrossSectionConv;
+
+	if (RadiationXRaySecondaryIon == TRUE) {
+	  nSecondaryHII = energies[bin] / 13.6;
+	  nSecondaryHeII = energies[bin] / 24.6;
+	}
+
+	/* Pre-calculate distances from cells to source */
+
+	for (dim = 0; dim < GridRank; dim++)
+	  for (i = 0, index = GridStartIndex[dim]; i < ActiveDims[dim]; 
+	       i++, index++) {
+
+	    // Calculate dr_i first, then square it
+	    ddr2[dim][i] = 
+	      fabs(CellLeftEdge[dim][index] + 0.5*CellWidth[dim][index] - 
+		   cstar->pos[dim]);
+	    ddr2[dim][i] = min(ddr2[dim][i], DomainWidth[dim]-ddr2[dim][i]);
+	    ddr2[dim][i] = ddr2[dim][i] * ddr2[dim][i];
+	  }
+
+	/* Loop over absorbers then cells */
+
+	double radius2_yz;
+
+	for (k = 0; k < ActiveDims[2]; k++) {
+	  for (j = 0; j < ActiveDims[1]; j++) {
+	    radius2_yz = ddr2[1][j] + ddr2[2][k];
+	    index = GRIDINDEX(0, j, k);
+	    for (i = 0; i < ActiveDims[0]; i++, index++) {
+	      r2_inv = 1.0 / (radius2_yz + ddr2[0][i]);
+	      if (RadiationXRaySecondaryIon == TRUE) {
+		xx = max(BaryonField[HIINum][index] / 
+			 (BaryonField[HINum][index] + BaryonField[HIINum][index]),
+			 1e-4);
+		heat_factor = 0.9971 * (1.0 - powf(1.0 - powf(xx, 0.2663f), 1.3163));
+		ion2_factor[0] = 0.3908 * nSecondaryHII * 
+		  powf(1 - powf(xx, 0.4092f), 1.7592f);
+		ion2_factor[1] = 0.0554 * nSecondaryHeII * 
+		  powf(1 - powf(xx, 0.4614f), 1.6660f);
+	      }
+	      for (a = 0; a < 3; a++) {
+		dkph = (float) (XrayLuminosity * CrossSections[a] / (4.0 * M_PI)) * r2_inv;
+		dE = (energies[bin] - EnergyThresholds[a]);
+		BaryonField[kphNum[a]][index] += dkph * ion2_factor[a];
+		BaryonField[gammaNum][index] += dkph * dE * heat_factor;
+	      } // END: absorbers (a)
+	    } // END: i-direction
+	  } // END: j-direction
+	} // END: k-direction
+
+      } // ENDFOR bin
+      
+    } // ENDFOR stars
+
+  } // ENDELSE ProblemType == 50
+
+  for (dim = 0; dim < GridRank; dim++)
+    delete [] ddr2[dim];
+
+  return SUCCESS;
+
+}

File src/enzo/Grid_AddXraysFromTree.C

+/***********************************************************************
+/
+/  ADD X-RAY EMISSION FROM STAR PARTICLES FROM A TREE
+/
+/  written by: John Wise
+/  date:       April, 2012
+/  modified1:
+/
+/  PURPOSE:
+/
+************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.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 "CosmologyParameters.h"
+#include "Star.h"
+
+#define MIN_OPENING_ANGLE 0.2  // 0.2 = arctan(11.3 deg)
+
+float CalculateLWFromTree(const FLOAT pos[], const float angle, 
+			  const SuperSourceEntry *Leaf, const float min_radius, 
+			  float result0);
+int FindSuperSourceByPosition(FLOAT *pos, SuperSourceEntry **result,
+			      int DEBUG);
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+
+int grid::AddXraysFromTree(void)
+{
+
+#ifdef UNUSED
+  const double pc = 3.086e18, clight = 3e10;
+
+  int i, j, k, index, dim, ci;
+  FLOAT pos[MAX_DIMENSION];
+  FLOAT radius2;
+  FLOAT innerFront, outerFront, innerFront2, outerFront2;
+  double Luminosity[MAX_ENERGY_BINS];
+  float energies[MAX_ENERGY_BINS], kdiss_r2;
+  double H2Luminosity, H2ISigma = 3.71e-18;
+
+  if (MyProcessorNumber != ProcessorNumber)
+    return SUCCESS;
+
+  /* Exit if there are no sources */
+
+  if (SourceClusteringTree == NULL)
+    return SUCCESS;
+
+  this->DebugCheck((char*) "Grid_AddH2Dissociation");
+
+  /* Get photo-ionization fields */
+
+  int kphHINum, kphHeINum, kphHeIINum, kdissH2INum;
+  int gammaNum;
+  IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, kphHeIINum, 
+				  kdissH2INum);
+
+  /* If using cosmology, get units. */
+
+  float TemperatureUnits, DensityUnits, LengthUnits, VelocityUnits, 
+    TimeUnits, aUnits = 1;
+
+  GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+	   &TimeUnits, &VelocityUnits, PhotonTime);
+
+  // Absorb the unit conversions into the cross-section
+  H2ISigma *= (double)TimeUnits / ((double)LengthUnits * (double)LengthUnits);
+
+  // Dilution factor (prevent breaking in the rate solver near the star)
+  float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+  float dilRadius2 = dilutionRadius * dilutionRadius;
+
+  // Convert from #/s to RT units
+  double LConv = (double) TimeUnits / pow(LengthUnits,3);
+  double LConv_inv = 1.0 / LConv;
+
+  /* Find sources in the tree that contribute to the cells */
+
+  SuperSourceEntry *Leaf;
+  float factor = LConv_inv * H2ISigma / (4.0 * M_PI);
+  float angle;
+
+  Leaf = SourceClusteringTree;
+
+  // We want to use the source separation instead of the merging
+  // radius.  The leaves store the merging radius (ClusteringRadius),
+  // so we multiply the angle by merge radius.
+  angle = MIN_OPENING_ANGLE * RadiativeTransferPhotonMergeRadius;
+
+  for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
+    pos[2] = CellLeftEdge[2][k] + 0.5*CellWidth[2][k];
+    for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+      pos[1] = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
+      index = GRIDINDEX_NOGHOST(GridStartIndex[0], j, k);
+      for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) {
+	pos[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
+
+	/* Find the leaves that have an opening angle smaller than
+	   the specified minimum and only include those in the
+	   calculation */
+
+	H2Luminosity = CalculateLWFromTree(pos, angle, Leaf, dilRadius2, 0);
+	BaryonField[kdissH2INum][index] = H2Luminosity * factor;
+
+      } // ENDFOR i
+    } // ENDFOR j
+  } // ENDFOR k
+#else
+  ENZO_FAIL("Optically thin X-rays calculation not implemented with tree acceleration.");
+#endif /* UNUSED */
+
+  return SUCCESS;
+
+}

File src/enzo/Grid_SetSubgridMarkerIsolatedBoundaries.C

+/***********************************************************************
+/
+/  GRID CLASS (MARK SUBGRID FOR ISOLATED BOUNDARIES)
+/
+/  written by: John Wise
+/  date:       April, 2012
+/  modified1:
+/
+/  PURPOSE:
+/
+/
+************************************************************************/
+
+#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"
+
+/* function prototypes */
+
+int grid::SetSubgridMarkerIsolatedBoundaries(void)
+{
+
+  /* Return if this grid is not on this processor. */
+
+  if (MyProcessorNumber != ProcessorNumber)
+    return SUCCESS;
+
+  /* declarations */
+    
+  int i, j, k, n, dim, field, index;
+  int SubgridStart[MAX_DIMENSION], SubgridEnd[MAX_DIMENSION];
+
+  /* Check if the outer cells lie outside the domain */
+
+  bool inside = true;
+  for (dim = 0; dim < MAX_DIMENSION; dim++)
+    inside &= CellLeftEdge[dim][0] > DomainLeftEdge[dim] &&
+      CellLeftEdge[dim][GridDimension[dim]-1] < DomainRightEdge[dim];
+  if (inside) return SUCCESS;
+
+  for (dim = 0; dim < MAX_DIMENSION; dim++) {
+    SubgridStart[dim] = 0;
+    SubgridEnd[dim] = 0;
+  }
+
+  /* Compute start and stop indices of the domain within this grid
+     (and check to make sure subgrid is within this grid). */
+
+  for (dim = 0; dim < GridRank; dim++) {
+
+    SubgridStart[dim] = nint(
+        (DomainLeftEdge[dim] - GridLeftEdge[dim])/CellWidth[dim][0]) 
+      + GridStartIndex[dim];
+    SubgridEnd[dim] = nint(
+	(DomainRightEdge[dim] - GridLeftEdge[dim])/CellWidth[dim][0]
+			       ) + GridStartIndex[dim] - 1;
+
+    SubgridStart[dim] = max(SubgridStart[dim], 0);
+    SubgridEnd[dim]   = min(SubgridEnd[dim], GridDimension[dim]-1);
+
+  }
+
+  /* Now that there is overlap mark this grid with a NULL pointer. */
+
+  int xdim, ydim, din;
+  for (dim = 0; dim < MAX_DIMENSION; dim++) {
+    xdim = (dim+1) % MAX_DIMENSION;
+    ydim = (dim+2) % MAX_DIMENSION;
+
+    // minus face
+    for (k = 0; k < SubgridStart[dim]; k++)
+      for (j = 0; j < GridDimension[ydim]; j++)
+	for (i = 0; i < GridDimension[xdim]; i++) {
+	  switch (dim) {
+	  case 0:
+	    index = GRIDINDEX_NOGHOST(k,i,j);
+	    break;
+	  case 1:
+	    index = GRIDINDEX_NOGHOST(j,k,i);
+	    break;
+	  case 2:
+	    index = GRIDINDEX_NOGHOST(i,j,k);
+	    break;
+	  }
+	  SubgridMarker[index] = NULL;
+	} // ENDFOR i
+    
+    // plus face
+    for (k = SubgridEnd[dim]+1; k < GridDimension[dim]; k++)
+      for (j = 0; j < GridDimension[ydim]; j++)
+	for (i = 0; i < GridDimension[xdim]; i++) {
+	  switch (dim) {
+	  case 0:
+	    index = GRIDINDEX_NOGHOST(k,i,j);
+	    break;
+	  case 1:
+	    index = GRIDINDEX_NOGHOST(j,k,i);
+	    break;
+	  case 2:
+	    index = GRIDINDEX_NOGHOST(i,j,k);
+	    break;
+	  }
+	  SubgridMarker[index] = NULL;
+	} // ENDFOR i
+
+  } // ENDFOR dim
+
+  return SUCCESS;
+  
+}

File src/enzo/Grid_Shine.C

 	  this_type = type_count;
 	  break;
 	}
-    }
+    } else if (RadiativeTransferOpticallyThinXray == TRUE) 
+      continue;
 
     if (DEBUG)
       printf("Shine: Photons/package[%"ISYM"]: %"GSYM" eV, %"GSYM", %"GSYM", %"GSYM", %"GSYM"\n", 

File src/enzo/Grid_WalkPhotonPackage.C

     g[dim] = GridStartIndex[dim] + 
       nint(floor((r[dim] - GridLeftEdge[dim]) / CellWidth[dim][0]));
     if (g[dim] < 0 || g[dim] >= GridDimension[dim]) {
-      //printf("Ray out of grid? g = %d %d %d\n", g[0], g[1], g[2]);
+      printf("Ray out of grid? g = %d %d %d\n", g[0], g[1], g[2]);
       DeleteMe = TRUE;
       return SUCCESS;
       //ENZO_FAIL("Ray out of grid?");
       return SUCCESS;
 
     // return in case we're out of photons
-    if ((*PP)->Photons < tiny_number || 
+    if ((*PP)->Photons < 1e-10*tiny_number || 
 	(*PP)->ColumnDensity > tau_delete) {
       if (DEBUG>1) {
 	fprintf(stderr, "PP-Photons: %"GSYM"  PP->Radius: %"GSYM

File src/enzo/Make.config.objects

 	Grid_AddH2Dissociation.o \
 	Grid_AddH2DissociationFromSources.o \
 	Grid_AddH2DissociationFromTree.o \
+	Grid_AddOpticallyThinXrays.o \
 	Grid_AddRadiationImpulse.o \
 	Grid_AddRadiationPressureAcceleration.o \
+	Grid_AddXraysFromSources.o \
+	Grid_AddXraysFromTree.o \
 	Grid_AllocateInterpolatedRadiation.o \
 	Grid_CommunicationSendPhotonPackages.o \
         Grid_ComputePhotonTimestep.o \
         Grid_SetSubgridMarkerFromParent.o \
         Grid_SetSubgridMarkerFromSibling.o \
         Grid_SetSubgridMarkerFromSubgrid.o \
+	Grid_SetSubgridMarkerIsolatedBoundaries.o \
         Grid_SubgridMarkerPostParallel.o \
         Grid_Shine.o \
         Grid_TransportPhotonPackages.o \

File src/enzo/PhotonGrid_Methods.h

    int SubgridMarkerPostParallel(HierarchyEntry **Grids[], int *NumberOfGrids);
    int SubgridMarkerPostParallelGZ(grid *Parent, HierarchyEntry **Grids[],
 				   int *NumberOfGrids);
+   int SetSubgridMarkerIsolatedBoundaries(void);
 
 /* Return Subgrid Marker for a position */
 
 
   int AddRadiationPressureAcceleration(void);
 
+/* Add optically thin X-ray field */
+
+  int AddXraysFromSources(Star *AllStars);
+  int AddXraysFromTree(void);
+  int AddOpticallyThinXrays(Star *AllStars, int NumberOfSources);
+
 /* Initialize ionized sphere around a source */
 
   int InitializeSource(RadiationSourceEntry *RS);

File src/enzo/RadiativeTransferParameters.h

 
 EXTERN int RadiativeTransferOpticallyThinH2;
 
+/* Flag to turn on a 1/r^2 X-ray radiation field */
+
+EXTERN int RadiativeTransferOpticallyThinXray;
+
 /* Periodic boundary conditions for the photon packages */
 
 EXTERN int RadiativeTransferPeriodicBoundary;

File src/enzo/RadiativeTransferReadParameters.C

   RadiativeTransferPropagationDistance        = 0.1;
   RadiativeTransferCoupledRateSolver          = TRUE;
   RadiativeTransferOpticallyThinH2            = TRUE;
+  RadiativeTransferOpticallyThinXray          = FALSE;
   RadiativeTransferFluxBackgroundLimit        = 0.01;
   RadiativeTransferSplitPhotonRadius          = FLOAT_UNDEFINED; // kpc
   RadiativeTransferRaysPerCell                = 5.1;
 		  &RadiativeTransferCoupledRateSolver);
     ret += sscanf(line, "RadiativeTransferOpticallyThinH2 = %"ISYM, 
 		  &RadiativeTransferOpticallyThinH2);
+    ret += sscanf(line, "RadiativeTransferOpticallyThinXray = %"ISYM, 
+		  &RadiativeTransferOpticallyThinXray);
     ret += sscanf(line, "RadiativeTransferPeriodicBoundary = %"ISYM, 
 		  &RadiativeTransferPeriodicBoundary);
     ret += sscanf(line, "RadiativeTransferSplitPhotonRadius = %"FSYM, 

File src/enzo/RadiativeTransferWriteParameters.C

 	  RadiativeTransferCoupledRateSolver);
   fprintf(fptr, "RadiativeTransferOpticallyThinH2          = %"ISYM"\n", 
 	  RadiativeTransferOpticallyThinH2);
+  fprintf(fptr, "RadiativeTransferOpticallyThinXray        = %"ISYM"\n", 
+	  RadiativeTransferOpticallyThinXray);
   fprintf(fptr, "RadiativeTransferFLDCallOnLevel           = %"ISYM"\n", 
 	  RadiativeTransferFLDCallOnLevel);
   fprintf(fptr, "RadiativeTransferPeriodicBoundary         = %"ISYM"\n", 

File src/enzo/RestartPhotons.C

   RadiativeTransferCoupledRateSolver = savedCoupledChemistrySolver;
   dtPhoton = SavedPhotonTimestep;
 
-  /* Optically thin Lyman-Werner (H2) radiation field */
+  /* Optically thin Lyman-Werner (H2) and X-ray radiation field */
 
   int NumberOfSources = 0;
-  Star *cstar = AllStars->NextStar;
-  while (cstar != NULL) {
-    cstar = cstar->NextStar;
-    NumberOfSources++;
+  if (AllStars != NULL) {
+    Star *cstar = AllStars->NextStar;
+    while (cstar != NULL) {
+      cstar = cstar->NextStar;
+      NumberOfSources++;
+    }
+  } else if (ProblemType == 50) {
+    RadiationSourceEntry *RS = GlobalRadiationSources->NextSource;
+    while (RS != NULL) {
+      RS = RS->NextSource;
+      NumberOfSources++;
+    }
   }
 
   if (RadiativeTransferOpticallyThinH2)
       for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel)
 	Temp->GridData->AddH2Dissociation(AllStars, NumberOfSources);
 
+  if (RadiativeTransferOpticallyThinXray)
+    for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
+      for (Temp = LevelArray[level]; Temp; Temp = Temp->NextGridThisLevel)
+	Temp->GridData->AddOpticallyThinXrays(AllStars, NumberOfSources);
+
   return SUCCESS;
 
 }