1. John Wise
  2. openmp

Commits

John Wise  committed 3c2f26a Merge

Merging.

  • Participants
  • Parent commits 6ad872f, 8ce7de7
  • Branches litemp

Comments (0)

Files changed (66)

File src/enzo/CommunicationPartitionGrid.C

View file
  
 // Function prototypes
  
-int CommunicationBroadcastValue(int *Value, int BroadcastProcessor);
 int Enzo_Dims_create(int nnodes, int ndims, int *dims);
 int LoadBalanceHilbertCurve(grid *GridPointers[], int NumberOfGrids, 
 			    int* &NewProcessorNumber);
+int CommunicationSyncNumberOfParticles(grid *GridPointer[], int NumberOfGrids);
 
 #define USE_OLD_CPU_DISTRIBUTION
 
   Unigrid = 0;
   if (debug) printf("Re-set Unigrid = 0\n");
  
-  /* Move Particles (while still on same processor). */
+  /* Move Particles (while still on same processor) and sync number of
+     particles. */
 
-  if (!ParallelRootGridIO)
-    if (OldGrid->MoveSubgridParticlesFast(gridcounter, SubGrids, TRUE) == FAIL) {
-      ENZO_FAIL("Error in grid->MoveSubgridParticlesFast.");
-    }
+  if (!ParallelRootGridIO) {
+    OldGrid->MoveSubgridParticlesFast(gridcounter, SubGrids, TRUE);
+    CommunicationSyncNumberOfParticles(SubGrids, gridcounter);
+  }
  
   int *PartitionProcessorNumbers = NULL;
   if (LoadBalancing == 4)
  
 	grid *NewGrid = ThisGrid->GridData;
  
-	/* Broadcast the number of particles to the other processors
-	   (OldGrid is assumed to be on the root processor). */
- 
-	if (NumberOfProcessors > 1) {
-
-	  int IntTemp = NewGrid->ReturnNumberOfParticles();
- 
-//          printf("NewGrid->ReturnNumberOfParticles: %"ISYM"\n", IntTemp);
- 
-	  CommunicationBroadcastValue(&IntTemp, ROOT_PROCESSOR);
-
-	  NewGrid->SetNumberOfParticles(IntTemp);
-
-//          printf("NG particle number set to %"ISYM"\n", IntTemp);
-
-	}
- 
 	/* Transfer from Old to New (which is still also on root processor) */
  
 	FLOAT Zero[] = {0,0,0};

File src/enzo/CommunicationSyncNumberOfParticles.C

View file
 {
 
   int i, idx;
-  //int *AllNumberOfParticles = new int[NumberOfGrids];
-  //int *AllNumberOfStars = new int[NumberOfGrids];
   int *buffer = new int[2*NumberOfGrids];
 
   for (i = 0, idx = 0; i < NumberOfGrids; i++, idx += 2)
 
 #ifdef USE_MPI
   CommunicationAllReduceValues(buffer, 2*NumberOfGrids, MPI_SUM);
-  //CommunicationAllReduceValues(AllNumberOfParticles, NumberOfGrids, MPI_SUM);
-  //CommunicationAllReduceValues(AllNumberOfStars, NumberOfGrids, MPI_SUM);
 #endif
 
   for (i = 0, idx = 0; i < NumberOfGrids; i++, idx += 2) {
   }
 
   delete [] buffer;
-  //delete [] AllNumberOfParticles;
-  //delete [] AllNumberOfStars;
 
   return SUCCESS;
 }
 
+/************************************************************************/
+
+int CommunicationSyncNumberOfParticles(grid *GridPointer[], int NumberOfGrids)
+{
+
+  int i, idx;
+  int *buffer = new int[2*NumberOfGrids];
+
+  for (i = 0, idx = 0; i < NumberOfGrids; i++, idx += 2)
+    if (GridPointer[i]->ReturnProcessorNumber() == MyProcessorNumber) {
+      buffer[idx] = GridPointer[i]->ReturnNumberOfParticles();
+      buffer[idx+1] = GridPointer[i]->ReturnNumberOfStars();
+    } else {
+      buffer[idx] = 0;
+      buffer[idx+1] = 0;
+    }
+
+#ifdef USE_MPI
+  CommunicationAllReduceValues(buffer, 2*NumberOfGrids, MPI_SUM);
+#endif
+
+  for (i = 0, idx = 0; i < NumberOfGrids; i++, idx += 2) {
+    GridPointer[i]->SetNumberOfParticles(buffer[idx]);
+    GridPointer[i]->SetNumberOfStars(buffer[idx+1]);
+  }
+
+  delete [] buffer;
+
+  return SUCCESS;
+}
+

File src/enzo/CreateSourceClusteringTree.C

View file
 //    if (SourceClusteringTree != NULL)
 //      DeleteSourceClusteringTree(SourceClusteringTree);
 
-  } // ENDIF SourceList == NULL (first time)  
+  } // ENDIF SourceList == NULL (first time)
 
   /* Calculate "center of light" first and assign it to the tree. */
 

File src/enzo/DetermineSubgridSizeExtrema.C

View file
   /* Now determine subgrid size parameters */
 
   int grids_per_proc = (level > MaximumStaticSubgridLevel) ?
-    OptimalSubgridsPerProcessor : 6;
+    OptimalSubgridsPerProcessor : 8;
 
   MaximumSubgridSize = NumberOfCells / 
     (NumberOfCores * grids_per_proc);

File src/enzo/EvolvePhotons.C

View file
 	  Grids[lvl][i]->GridData->FinalizeRadiationFields();
     END_PERF(8);
 
-    /* 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)
 #pragma omp parallel for schedule(guided)
 	for (i = 0; i < nGrids[lvl]; i++)
 	  Grids[lvl][i]->GridData->AddH2Dissociation(AllStars, NumberOfSources);
+    if (RadiativeTransferOpticallyThinXray)
+      for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
+#pragma omp parallel for schedule(guided)
+	for (i = 0; i < nGrids[lvl]; i++)
+	  Grids[lvl][i]->GridData->
+	    AddOpticallyThinXrays(AllStars, NumberOfSources);
     END_PERF(10);
 
     START_PERF();
 #pragma omp parallel for schedule(guided)
 	for (i = 0; i < nGrids[lvl]; i++)
 	  if (Grids[lvl][i]->GridData->RadiationPresent() == TRUE) {
-
 	    Grids[lvl][i]->GridData->SolveRateAndCoolEquations
 	      (RTCoupledSolverIntermediateStep);
-
 	  } /* ENDIF radiation */
     END_PERF(9);
 

File src/enzo/FindSuperSourceByPosition.C

View file
   __m128 x = _mm_set_ss(*__x);
   __m128 recip = _mm_rsqrt_ss(x);
   _mm_store_ss(__outrsqrt, recip);
-// _m128* precip = (__m128 *)__outrsqrt;
-// *precip = _mm_mul_ss(_mm_set_ss(0.5f), _mm_add_ss(recip, _mm_rcp_ss(_mm_mul_ss(x, recip))));
 }
 #endif
 

File src/enzo/Grid.h

View file
 
    int FlagCellsToBeRefinedByMetallicity(int level);
 
+/* Flag all cells which have more than a specified metal mass */
+
+   int FlagCellsToBeRefinedByMetalMass(int level);
 
 /* Flagging all cell adjacent to a previous flagged cell.  Also, remove all
    Flagged cells in the boundary zones and within one zone of the boundary. */

File src/enzo/Grid_AddOpticallyThinXrays.C

View file
+/***********************************************************************
+/
+/  (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

View file
+/***********************************************************************
+/
+/  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

View file
+/***********************************************************************
+/
+/  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_ComputeCoolingTime.C

View file
         int *idual, int *ispecies, int *imetal, int *imcool, int *idust, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke, int *ih2co,
 	int *ipiht, int *igammah,
-	float *dt, float *aye, float *temstart, float *temend,
+	float *dt, float *aye, float *redshift, float *temstart, float *temend,
 	float *utem, float *uxyz, float *uaye, float *urho, float *utim,
 	float *eta1, float *eta2, float *gamma, float *z_solar,
 	float *ceHIa, float *ceHeIa, float *ceHeIIa, float *ciHIa, float *ciHeIa,
        &GridRank, GridStartIndex, GridStartIndex+1, GridStartIndex+2,
        GridEndIndex, GridEndIndex+1, GridEndIndex+2,
        &CoolData.ih2co, &CoolData.ipiht, &PhotoelectricHeating,
-       &dtFixed, &afloat, &CoolData.TemperatureStart,
-       &CoolData.TemperatureEnd,
+       &dtFixed, &afloat, &RadiationFieldRedshift,
+       &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
        &TemperatureUnits, &LengthUnits, &aUnits, &DensityUnits, &TimeUnits,
        &DualEnergyFormalismEta1, &DualEnergyFormalismEta2, &Gamma,
        &CoolData.SolarMetalFractionByMass,

File src/enzo/Grid_ComputeTimeStep.C

View file
     ENZO_FAIL("Error in GetUnits.\n");
   }
 
+  float dtStar = huge_number;
+#ifdef UNUSED
   float mindtNOstars;  // Myr
   const int NumberOfStepsInLifetime = 5;
-  float dtStar = huge_number;
 
   if (STARFEED_METHOD(POP3_STAR))
     mindtNOstars = 3;  // Myr
       dtStar = 3.1557e13*mindtNOstars/TimeUnits;
 
   dt = min(dt, dtStar);
+#endif /* UNUSED */
 
 
 
       printf("Cond = %"ESYM" ",(dtConduction));
     if (UseGasDrag)
       printf("Drag = %"ESYM" ",(dtGasDrag));
+#ifdef TRANSFER
+    if (dtStar < huge_number)
+      printf("Stars = %"FSYM" ", dtStar);
+    if (RadiationPressure)
+      printf("RP = %"FSYM" ", dtRadPressure);
+    if (TimestepSafetyVelocity > 0)
+      printf("Saf = %"FSYM" ", dtSafetyVelocity);
+#endif /* TRANSFER */      
     printf(")\n");
   }
  

File src/enzo/Grid_CopyZonesFromGrid.C

View file
  
   /* declarations */
  
-  int i,j,k,field,dim;
+  int i, j, k, dim, loop, iLoop, field, otherindexB;
+  float vx, vy, vz, v2, rho, bx, by, bz, b2;
 
   bool shiftPos, shiftNeg; float delta; FLOAT L;
 
   if (ShearingBoundaryDirection!=-1){
     L=(DomainRightEdge[ShearingBoundaryDirection]-DomainLeftEdge[ShearingBoundaryDirection]);
 
-   
-    bool noMove=false;
-
     delta=L*AngularVelocity*VelocityGradient;
   
     //printf("L: %"GSYM" Delta: %"GSYM" %"GSYM" (%"GSYM" %"GSYM")\n", L, delta, delta, AngularVelocity, VelocityGradient);
 
 
   //Shearing Boundary Variables
-  float rho, vx, vy, vz, v2, b2, bx, by, bz=0.0;
   int thisindex, otherindex=0;
   FLOAT a,b;  FLOAT val1=-9999;FLOAT val2=-9999;
   
 
   int addDim[3] = {1, OtherDim[0], OtherDim[0]*OtherDim[1]};
   int velocityTypes[3]={Velocity1, Velocity2, Velocity3};
-  int Zero[3] = {0,0,0};
 
   if (!isShearing)
-    //#pragma omp parallel for schedule(static)
+//#pragma omp parallel for schedule(static)
     for (field = 0; field < NumberOfBaryonFields; field++)
       FORTRAN_NAME(copy3drel)(OtherGrid->BaryonField[field], BaryonField[field],
 			      Dim, Dim+1, Dim+2,
 	    (k + StartOther[2])*OtherDim[0]*OtherDim[1];
 	  for (i = 0; i < Dim[0]; i++, thisindex++, otherindex++){
 
-	    int otherindexB=otherindex+ addDim[ShearingVelocityDirection];
+	    otherindexB=otherindex+ addDim[ShearingVelocityDirection];
 	    
 	    val1=OtherGrid->BaryonField[field][otherindex];
 	    val2=OtherGrid->BaryonField[field][otherindexB];
 	    
 	    if (DualEnergyFormalism==0 && FieldType[field]==TotalEnergy) {
-	      for (int loop=0; loop<=1; loop++){
-		int iLoop=otherindex;
+	      for (loop=0; loop < 2; loop++){
+		iLoop=otherindex;
 		if (loop==1) iLoop= otherindexB;
 
-		float vx, vy, vz, v2, rho;
 		v2=0.0;
 		rho= OtherGrid->BaryonField[iden][iLoop];
 		vx= OtherGrid->BaryonField[ivx][iLoop];
 		else vz=0.0;
 		v2=vx*vx+vy*vy+vz*vz;
 		
-		float bx, by, bz, b2;
 		b2=0.0;
 		if (useMHD) {
 		  bx= OtherGrid->BaryonField[iBx][iLoop];
 	thisindex = (0 + Start[0]) + (j + Start[1])*GridDimension[0] +
 	  (k + Start[2])*GridDimension[0]*GridDimension[1];
 	for (i = 0; i < Dim[0]; i++, thisindex++){
-	  float vx, vy, vz, v2, rho, bx, by, bz, b2;
 	  rho= BaryonField[iden][thisindex];  
 	  vx= BaryonField[ivx][thisindex];
 	  vy= BaryonField[ivy][thisindex];  

File src/enzo/Grid_FlagCellsToBeRefinedByMetalMass.C

View file
+/***********************************************************************
+/
+/  GRID CLASS (FLAG CELLS TO BE REFINED BY METAL MASS)
+/
+/  written by: John Wise
+/  date:       March, 2012
+/  modified1:  
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    number of flagged cells, or -1 on failure
+/
+************************************************************************/
+ 
+#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"
+ 
+int grid::FlagCellsToBeRefinedByMetalMass(int level)
+{
+
+  int i, dim, size, method;
+  const float Zsun = 0.0204;
+ 
+  /* Return if this grid is not on this processor. */
+ 
+  if (MyProcessorNumber != ProcessorNumber)
+    return SUCCESS;
+ 
+  /* error check */
+
+  if (FlaggingField == NULL) {
+    fprintf(stderr, "Flagging Field is undefined.\n");
+    return -1;
+  }
+
+  /* Search for mass flagging field */
+
+  int MassFlaggingMethod = INT_UNDEFINED;
+  for (method = 0; method < MAX_FLAGGING_METHODS; method++)
+    if (CellFlaggingMethod[method] == 2) {
+      MassFlaggingMethod = method;
+      break;
+    }
+
+  if (MassFlaggingMethod == INT_UNDEFINED)
+    ENZO_FAIL("Metal mass refinement requires baryon refinement.");
+
+  /* compute size */
+
+  size = 1;
+  for (dim = 0; dim < GridRank; dim++)
+    size *= GridDimension[dim];
+
+  /* Compute cell volume */
+ 
+  float CellVolume = 1.0;
+  for (dim = 0; dim < GridRank; dim++)
+    CellVolume *= CellWidth[dim][0];
+ 
+  /* Compute the ModifiedMinimumMass */
+ 
+  float ModifiedMinimumMassForRefinement =
+    (MetallicityForRefinement*Zsun) * 
+    MinimumMassForRefinement[MassFlaggingMethod] * 
+    POW(RefineBy, level*MinimumMassForRefinementLevelExponent[MassFlaggingMethod]);
+
+  /* Metal cooling field numbers. */
+
+  int MetalNum = 0, SNColourNum = 0;
+  int MetalFieldPresent = FALSE;
+
+  // First see if there's a metal field (so we can conserve species in
+  // the solver)
+  MetalNum = FindField(Metallicity, FieldType, NumberOfBaryonFields);
+  SNColourNum = FindField(SNColour, FieldType, NumberOfBaryonFields);
+  MetalFieldPresent = (MetalNum != -1 || SNColourNum != -1);
+
+  // Double check if there's a metal field when we have metal cooling
+  if (MetalCooling && MetalFieldPresent == FALSE)
+    ENZO_FAIL("ERROR: No metal field found. "
+	      "Restart and turn OFF MetalMass refinement.");
+
+  /* If both metal fields (Pop I/II and III) exist, create a field
+     that contains their sum */
+
+  float *MetalPointer;
+  float *TotalMetals = NULL;
+
+  if (MetalNum != -1 && SNColourNum != -1) {
+    TotalMetals = new float[size];
+    for (i = 0; i < size; i++)
+      TotalMetals[i] = BaryonField[MetalNum][i] + BaryonField[SNColourNum][i];
+    MetalPointer = TotalMetals;
+  } // ENDIF both metal types
+  else {
+    if (MetalNum != -1)
+      MetalPointer = BaryonField[MetalNum];
+    else if (SNColourNum != -1)
+      MetalPointer = BaryonField[SNColourNum];
+  } // ENDELSE both metal types
+
+  /* Flag points */
+
+  for (i = 0; i < size; i++)
+    FlaggingField[i] += 
+      (CellVolume * MetalPointer[i] > ModifiedMinimumMassForRefinement) ? 1 : 0;
+ 
+  /* Count number of flagged Cells. */
+ 
+  int NumberOfFlaggedCells = 0;
+  for (i = 0; i < size; i++)
+    if (FlaggingField[i] > 0)
+      NumberOfFlaggedCells++;
+ 
+  delete [] TotalMetals;
+ 
+  return NumberOfFlaggedCells;
+ 
+}

File src/enzo/Grid_Group_WriteGrid.C

View file
     /* If requested, compute and output the dust temperature field 
        as well since its such a pain to compute after the fact. */
 
-    if (OutputDustTemperature) {
+    if (OutputDustTemperature != FALSE) {
  
       /* Get temperature field if we do not already have it. */
 

File src/enzo/Grid_MergePausedPhotonPackages.C

File contents unchanged.

File src/enzo/Grid_PhotonTestInitializeGrid.C

View file
 			     float PhotonTestInitialFractionH2II,
 			     int RefineByOpticalDepth,
 			     int TotalRefinement,
-			     char *DensityFilename)
+			     char *DensityFilename,
+			     char *HIIFractionFilename,
+			     char *HeIIFractionFilename,
+			     char *HeIIIFractionFilename,
+			     char *TemperatureFilename)
 {
   /* declarations */
 
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
     DINum, DIINum, HDINum,  kphHINum, gammaNum, kphHeINum,
     kphHeIINum, kdissH2INum, RPresNum1, RPresNum2, RPresNum3; 
-  float *density_field = NULL;
+  float *density_field = NULL, *HII_field = NULL, *HeII_field = NULL, 
+    *HeIII_field = NULL, *Temperature_field = NULL;
 
   /* create fields */
   NumberOfBaryonFields = 0;
   const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
                pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16;
   float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, 
-    VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
+    VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6, mu_data;
 
   FLOAT a, dadt, ExpansionFactor = 1;
   GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
     dataset_name = strtok(NULL, delim);
 
     file_id = H5Fopen(data_filename, H5F_ACC_RDONLY, H5P_DEFAULT);
-    if (file_id == -1) ENZO_FAIL("Error closing density field.");
+    if (file_id == -1) ENZO_FAIL("Error opening density field.");
     for (dim = 0; dim < MAX_DIMENSION; dim++)
       OutDims[GridRank-dim-1] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
     density_field = new float[active_size];
 		       HDF5_REAL, density_field, FALSE, NULL, NULL);
     h5error = H5Fclose(file_id);
     if (h5error == -1) ENZO_FAIL("Error closing density field.");
+
+    /* HII, HeI, HeII, HeIII density fields */
+
+    if (HIIFractionFilename != NULL) {
+      data_filename = strtok(HIIFractionFilename, delim);
+      dataset_name = strtok(NULL, delim);
+      file_id = H5Fopen(data_filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+      if (file_id == -1) ENZO_FAIL("Error opening HII field.");
+      for (dim = 0; dim < MAX_DIMENSION; dim++)
+	OutDims[GridRank-dim-1] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
+      HII_field = new float[active_size];
+      this->read_dataset(GridRank, OutDims, dataset_name, file_id,
+			 HDF5_REAL, HII_field, FALSE, NULL, NULL);
+      h5error = H5Fclose(file_id);
+      if (h5error == -1) ENZO_FAIL("Error closing HII field.");
+    }
+
+    if (HeIIFractionFilename != NULL) {
+      data_filename = strtok(HeIIFractionFilename, delim);
+      dataset_name = strtok(NULL, delim);
+      file_id = H5Fopen(data_filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+      if (file_id == -1) ENZO_FAIL("Error opening HeII field.");
+      for (dim = 0; dim < MAX_DIMENSION; dim++)
+	OutDims[GridRank-dim-1] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
+      HeII_field = new float[active_size];
+      this->read_dataset(GridRank, OutDims, dataset_name, file_id,
+			 HDF5_REAL, HeII_field, FALSE, NULL, NULL);
+      h5error = H5Fclose(file_id);
+      if (h5error == -1) ENZO_FAIL("Error closing HeII field.");
+    }
+
+    if (HeIIIFractionFilename != NULL) {
+      data_filename = strtok(HeIIIFractionFilename, delim);
+      dataset_name = strtok(NULL, delim);
+      file_id = H5Fopen(data_filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+      if (file_id == -1) ENZO_FAIL("Error opening HeIII field.");
+      for (dim = 0; dim < MAX_DIMENSION; dim++)
+	OutDims[GridRank-dim-1] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
+      HeIII_field = new float[active_size];
+      this->read_dataset(GridRank, OutDims, dataset_name, file_id,
+			 HDF5_REAL, HeIII_field, FALSE, NULL, NULL);
+      h5error = H5Fclose(file_id);
+      if (h5error == -1) ENZO_FAIL("Error closing HeIII field.");
+    }
+
+    if (TemperatureFilename != NULL) {
+      data_filename = strtok(TemperatureFilename, delim);
+      dataset_name = strtok(NULL, delim);
+      file_id = H5Fopen(data_filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+      if (file_id == -1) ENZO_FAIL("Error opening temperature field.");
+      for (dim = 0; dim < MAX_DIMENSION; dim++)
+	OutDims[GridRank-dim-1] = GridEndIndex[dim] - GridStartIndex[dim] + 1;
+      Temperature_field = new float[active_size];
+      this->read_dataset(GridRank, OutDims, dataset_name, file_id,
+			 HDF5_REAL, Temperature_field, FALSE, NULL, NULL);
+      h5error = H5Fclose(file_id);
+      if (h5error == -1) ENZO_FAIL("Error closing temperature field.");
+    }
+
+
   } // ENDIF DensityFilename
 
   /* Loop over the mesh. */
 
 	/* Loop over spheres. */
 
-	if (density_field != NULL && 
-	    i >= GridStartIndex[0] && i <= GridEndIndex[0] &&
+	if (i >= GridStartIndex[0] && i <= GridEndIndex[0] &&
 	    j >= GridStartIndex[1] && j <= GridEndIndex[1] &&
 	    k >= GridStartIndex[2] && k <= GridEndIndex[2]) {
 	  cindex = (i-GridStartIndex[0]) + ActiveDims[0] *
 	    ((j-GridStartIndex[1]) + (k-GridStartIndex[2])*ActiveDims[1]);
-	  density = density_field[cindex];
-	} else {
-	  density = 1.0;
+	  if (density_field != NULL)
+	    density = max(density_field[cindex], 1e-6);
+	  else
+	    density = 1.0;
+	  if (HII_field != NULL)
+	    HII_Fraction = HII_field[cindex];
+	  else
+	    HII_Fraction = PhotonTestInitialFractionHII;
+	  if (HeII_field != NULL)
+	    HeII_Fraction = HeII_field[cindex];
+	  else
+	    HeII_Fraction = PhotonTestInitialFractionHeII;
+	  if (HeIII_field != NULL)
+	    HeIII_Fraction = HeIII_field[cindex];
+	  else
+	    HeIII_Fraction = PhotonTestInitialFractionHeIII;
+	  if (Temperature_field != NULL)
+	    temperature = temp1 = Temperature_field[cindex];
+	  else
+	    temperature = temp1 = InitialTemperature;
 	}
 
-	temperature = temp1 = InitialTemperature;
+	H2I_Fraction = PhotonTestInitialFractionH2I;
 	sigma = sigma1 = 0;
 	colour = 1.0e-10;
 	for (dim = 0; dim < MAX_DIMENSION; dim++)
 	  Velocity[dim] = 0;
-	HII_Fraction = PhotonTestInitialFractionHII;
-	HeII_Fraction = PhotonTestInitialFractionHeII;
-	HeIII_Fraction = PhotonTestInitialFractionHeIII;
-	H2I_Fraction = PhotonTestInitialFractionH2I;
 
 	for (sphere = 0; sphere < NumberOfSpheres; sphere++) {
 
 	  BaryonField[ivel+dim][n] = Velocity[dim] + UniformVelocity[dim];
 	
 	/* Set energy (thermal and then total if necessary). */
-	
+
+	if (MultiSpecies) {
+	  mu_data =  
+	    0.25*(BaryonField[HeINum][n]  + BaryonField[HeIINum][n] +
+		  BaryonField[HeIIINum][n]                        ) +
+	    BaryonField[HINum][n]   + BaryonField[HIINum][n]  +
+	    BaryonField[DeNum][n];
+	  if (MultiSpecies > 1)
+	    mu_data += BaryonField[HMNum][i]   +
+	      0.5*(BaryonField[H2INum][i]  + BaryonField[H2IINum][i]);
+	  mu_data = BaryonField[0][n] / mu_data;
+	} else
+	  mu_data = mu;
+
 	BaryonField[1][n] = temperature/TemperatureUnits/
-	  ((Gamma-1.0)*mu);
+	  ((Gamma-1.0)*mu_data);
 	
 	if (DualEnergyFormalism)
 	  BaryonField[2][n] = BaryonField[1][n];
     printf("PhotonTestInitialize: DM NumberOfParticles = %"ISYM"\n", 
 	   NumberOfParticles);
 
-  if (density_field != NULL)
-    delete [] density_field;
+  delete [] density_field;
+  delete [] HII_field;
+  delete [] HeII_field;
+  delete [] HeIII_field;
+  delete [] Temperature_field;
   
   return SUCCESS;
 }

File src/enzo/Grid_ProjectToPlane2.C

View file
 
     first_field = temperature;
     second_field = BaryonField[0];
-    ProjType = 3;
-    ConversionFactor = DensityConversion*DensityConversion;
+    ProjType = 4;
+    ConversionFactor = DensityConversion;//*DensityConversion;
 
     FORTRAN_NAME(projplane)(first_field, second_field,
 			    BaryonField[NumberOfBaryonFields], &One,
 
   if (NumberOfBaryonFields > 0) {
 
-    for (i = 0; i < size; i++)
-      temp_field[i] = BaryonField[0][i] * BaryonField[0][i];
+//    for (i = 0; i < size; i++)
+//      temp_field[i] = BaryonField[0][i] * BaryonField[0][i];
     
-    first_field = temp_field;
-    second_field = NULL;
-    ProjType = 1;
-    ConversionFactor = DensityConversion*DensityConversion;
+    first_field = BaryonField[0];
+    second_field = BaryonField[0];
+    ProjType = 4;
+    ConversionFactor = DensityConversion;//*DensityConversion;
 
     FORTRAN_NAME(projplane)(first_field, second_field,
                              BaryonField[NumberOfBaryonFields], &One,
     first_field = temp_field;
     second_field = BaryonField[0];
     //second_field = NULL;
-    ProjType = 3;
-    ConversionFactor = DensityConversion*DensityConversion/CoolData.SolarMetalFractionByMass;
+    ProjType = 4;
+    ConversionFactor = DensityConversion/CoolData.SolarMetalFractionByMass;
     //ConversionFactor = One;
 
     FORTRAN_NAME(projplane)(first_field, second_field,
     
     first_field = temp_field;
     second_field = BaryonField[0];
-    ProjType = 3;
-    ConversionFactor = DensityConversion*DensityConversion;
+    ProjType = 4;
+    ConversionFactor = DensityConversion;//*DensityConversion;
 
     FORTRAN_NAME(projplane)(first_field, second_field,
                              BaryonField[NumberOfBaryonFields], &One,
 
     first_field = temp_field;
     second_field = BaryonField[0];
-    ProjType = 3;
-    ConversionFactor = DensityConversion*DensityConversion;
+    ProjType = 4;
+    ConversionFactor = DensityConversion;//*DensityConversion;
 
     FORTRAN_NAME(projplane)(first_field, second_field,
                              BaryonField[NumberOfBaryonFields], &One,

File src/enzo/Grid_SetFlaggingField.C

View file
        is automatically turned if method #8 is specified. */
 
     break;
- 
+
+    /* ==== METHOD 19: Refine on metal mass ==== */
+
+  case 19:
+    NumberOfFlaggedCells = this->FlagCellsToBeRefinedByMetalMass(level);
+    if (NumberOfFlaggedCells < 0) {
+      fprintf(stderr, "Error in grid->FlagCellsToBeRefinedByMetalMass.\n");
+      return FAIL;
+    }
+    break;
  
     /* ==== METHOD 100: UNDO REFINEMENT IN SOME REGIONS ==== */
  
  
   }
 
+  if (debug1 && NumberOfFlaggedCells > 0)
+    printf("SetFlaggingField[method = %"ISYM"]: NumberOfFlaggedCells = %"ISYM".\n",
+	   CellFlaggingMethod[method], NumberOfFlaggedCells);
+
   } // ENDFOR methods
  
   /* End of Cell flagging criterion routine                              */
   counter[4]++;
   timer[4] += NumberOfFlaggedCells;
 #endif /* MPI_INSTRUMENTATION */
- 
-  if (debug1)
-
-    printf("SetFlaggingField[method = %"ISYM"]: NumberOfFlaggedCells = %"ISYM".\n",
-	   method, NumberOfFlaggedCells);
- 
+  
   return SUCCESS;
  
 }

File src/enzo/Grid_SetSubgridMarkerIsolatedBoundaries.C

View file
+/***********************************************************************
+/
+/  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

View file
 	  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_SolveRadiativeCooling.C

View file
         int *idual, int *ispecies, int *imetal, int *imcool, int *idust, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke, int *imax, int *ih2co,
 	int *ipiht,
-	float *dt, float *aye, float *temstart, float *temend,
+	float *dt, float *aye, float *redshift, float *temstart, float *temend,
 	float *utem, float *uxyz, float *uaye, float *urho, float *utim,
 	float *eta1, float *eta2, float *gamma, float *z_solar,
 	float *ceHIa, float *ceHeIa, float *ceHeIIa, float *ciHIa,
        &GridRank, GridStartIndex, GridStartIndex+1, GridStartIndex+2,
        GridEndIndex, GridEndIndex+1, GridEndIndex+2, &MaxDimension,
        &CoolData.ih2co, &CoolData.ipiht,
-       &dtFixed, &afloat, &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
+       &dtFixed, &afloat, &RadiationFieldRedshift,
+       &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
        &TemperatureUnits, &LengthUnits, &aUnits, &DensityUnits, &TimeUnits,
        &DualEnergyFormalismEta1, &DualEnergyFormalismEta2, &Gamma, 
        &CoolData.SolarMetalFractionByMass,

File src/enzo/Grid_SolveRateAndCoolEquations.C

View file
         int *idual, int *ispecies, int *imetal, int *imcool, int *idust, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke, int *imax, int *ih2co, 
 	int *ipiht, int *igammah,
-	float *dt, float *aye, float *temstart, float *temend,
+	float *dt, float *aye, float *redshift, float *temstart, float *temend,
 	float *utem, float *uxyz, float *uaye, float *urho, float *utim,
 	float *eta1, float *eta2, float *gamma, float *fh, float *dtoh,
 	float *z_solar,
     &GridRank, GridStartIndex, GridStartIndex+1, GridStartIndex+2, 
     GridEndIndex, GridEndIndex+1, GridEndIndex+2, &MaxDimension,
     &CoolData.ih2co, &CoolData.ipiht, &PhotoelectricHeating,
-    &dtCool, &afloat, &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
+    &dtCool, &afloat, &RadiationFieldRedshift, 
+    &CoolData.TemperatureStart, &CoolData.TemperatureEnd,
     &TemperatureUnits, &LengthUnits, &aUnits, &DensityUnits, &TimeUnits,
     &DualEnergyFormalismEta1, &DualEnergyFormalismEta2, &Gamma,
     &CoolData.HydrogenFractionByMass, &CoolData.DeuteriumToHydrogenRatio,

File src/enzo/Grid_TransferSubgridParticles.C

View file
 	   block, so we know where to start inserting them in the move
 	   list. */
 
+#ifdef _OPENMP
 	move_per_thread = nmove / CoresPerProcess;
 	for (i = 0; i < CoresPerProcess+1; i++)
 	  ThreadCount[i] = 0;
 	      n1++;
 	    } // ENDIF move to subgrid
 	  } // ENDFOR particles
-
+#endif /* _OPENMP */
       } // ENDELSE SingleThread
 
     } // END parallel

File src/enzo/Grid_WalkPhotonPackage.C

View file
       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/Grid_WriteGrid.C

View file
     /* If requested compute and output the dust temperature field
        as well since its such a pain to compute after the fact. */
   
-    if (OutputDustTemperature) {
+    if (OutputDustTemperature != FALSE) {
  
       /* Get temperature field if we do not already have it. */
 

File src/enzo/InterpretCommandLine.C

View file
 			 int &WritePotentialOnly,
 			 int &SmoothedDarkMatterOnly,
 			 int &WriteCoolingTimeOnly,
+			 int &WriteDustTemperatureOnly,
 			 int MyProcessorNumber)
 {
  
 	WriteCoolingTimeOnly = TRUE;
 	break;
 
+	/* Add dust temperature to data */
+
+      case 'D':
+	WriteDustTemperatureOnly = TRUE;
+	break;
+
 	/* debug */
  
       case 'd':
 	          "      -M (Write smoothed DM field only)\n"
 	          "      -F(riends-of-friends halo finder only)\n"
 	          "      -C(ooling time write only)\n"
+	          "      -D(ust temperature write only)\n"
                   "      -h(elp)\n"
 	          "      -i(nformation output)\n"
 	          "      -V (show compiler options and flags)\n"

File src/enzo/Make.config.objects

View file
 	Grid_FlagCellsToBeRefinedByTotalJeansLength.o \
 	Grid_FlagCellsToBeRefinedByMass.o \
 	Grid_FlagCellsToBeRefinedByMetallicity.o \
+	Grid_FlagCellsToBeRefinedByMetalMass.o \
 	Grid_FlagCellsToBeRefinedByMustRefineParticles.o \
 	Grid_FlagCellsToBeRefinedByMustRefineRegion.o \
 	Grid_FlagCellsToBeRefinedByResistiveLength.o \
 	OneZoneFreefallTestInitialize.o \
         OutputAsParticleData.o \
 	OutputCoolingTimeOnly.o \
+	OutputDustTemperatureOnly.o \
         OutputFromEvolveLevel.o\
         OutputLevelInformation.o \
 	OutputPotentialFieldOnly.o \
 	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/NestedCosmologySimulationInitialize.C

View file
     // Go down to the grid(s) on the next level
  
     CurrentGrid = CurrentGrid->NextGridNextLevel;
- 
+
+    CommunicationBarrier();
+     
   } // end loop over initial grid levels
 
   // Create tracer particles (on top grid)

File src/enzo/New_Grid_WriteGrid.C

View file
     /* If requested, compute and output the dust temperature field 
        as well since its such a pain to compute after the fact. */
  
-    if (OutputDustTemperature) {
+    if (OutputDustTemperature != FALSE) {
  
       /* Get temperature field if we do not already have it. */
 

File src/enzo/OutputDustTemperatureOnly.C

View file
+/***********************************************************************
+/
+/  COMPUTE THE DUST TEMPERATURE AND OUTPUT
+/
+/  written by: John Wise
+/  date:       October, 2011
+/  modified1:  
+/
+/  PURPOSE:
+/
+************************************************************************/
+#include "preincludes.h"
+
+#include <stdlib.h>
+#include <string.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"
+#include "Hierarchy.h"
+#include "TopGridData.h"
+#ifdef TRANSFER
+#include "ImplicitProblemABC.h"
+#endif
+#include "CosmologyParameters.h"
+
+int Group_WriteAllData(char *basename, int filenumber,
+		       HierarchyEntry *TopGrid, TopGridData &MetaData,
+		       ExternalBoundary *Exterior, 
+#ifdef TRANSFER
+		       ImplicitProblemABC *ImplicitSolver,
+#endif
+		       FLOAT WriteTime = -1, 
+		       int CheckpointDump = FALSE);
+
+int RebuildHierarchy(TopGridData *MetaData,
+		     LevelHierarchyEntry *LevelArray[], int level);
+#ifdef TRANSFER
+int RadiativeTransferReadParameters(FILE *fptr);
+#endif
+
+int OutputDustTemperatureOnly(char *ParameterFile,
+			      LevelHierarchyEntry *LevelArray[], 
+			      HierarchyEntry *TopGrid,
+			      TopGridData &MetaData,
+			      ExternalBoundary &Exterior
+#ifdef TRANSFER
+			      , ImplicitProblemABC *ImplicitSolver
+#endif
+			  )
+{
+
+  int i, level;
+  const float When = 0.5;
+  LevelHierarchyEntry *Temp;
+  HierarchyEntry **Grids;
+
+  /* Exit if OutputDustTemperature is already on. */
+
+  if (OutputDustTemperature == TRUE) {
+    if (MyProcessorNumber == ROOT_PROCESSOR)
+      printf("OutputDustTemperature is already on.  "
+	     "Change to 0 in %s to recompute.\n", ParameterFile);
+    return SUCCESS;
+  }
+
+  if (MetalCooling == FALSE)
+    if (MyProcessorNumber == ROOT_PROCESSOR)
+      printf("Metal cooling is OFF.  Dust temperatures may be inaccurate.\n");
+
+  /* Determine the parameter name prefix */
+
+  char *cptr;
+  int DumpNumber;
+  char DumpName[MAX_LINE_LENGTH];
+  char NumberString[MAX_LINE_LENGTH];
+  if ( (cptr = strstr(ParameterFile, MetaData.DataDumpName)) ) {
+    strcpy(DumpName, MetaData.DataDumpName);
+  }
+  else if ( (cptr = strstr(ParameterFile, MetaData.RestartDumpName)) ) {
+    strcpy(DumpName, MetaData.RestartDumpName);
+  }
+  else if ( (cptr = strstr(ParameterFile, MetaData.RedshiftDumpName)) ) {
+    strcpy(DumpName, MetaData.RedshiftDumpName);
+  }
+  else
+    ENZO_FAIL("Cannot determine output type.");
+
+  /* Extract output number */
+  
+  strncpy(NumberString, ParameterFile+2, 4);
+  NumberString[4] = '\0';
+  DumpNumber = atoi(NumberString);
+
+  /* If we're not using parallel root grid I/O and we're parallel, we
+     need to rebuild the hierarchy with the multiple root grids. */
+
+  if (!ParallelRootGridIO && NumberOfProcessors > 1)
+    RebuildHierarchy(&MetaData, LevelArray, 0);
+
+  /* Initialize radiative transfer parameters, if needed */
+
+#ifdef TRANSFER
+  FILE *fptr;
+
+  if ((fptr = fopen(ParameterFile, "r")) == NULL) {
+    ENZO_VFAIL("Error opening ParameterFile %s\n", ParameterFile)
+  }
+
+  RadiativeTransferReadParameters(fptr);
+
+  fclose(fptr);
+#endif /* TRANSFER */
+
+  // Negative number to indicate that this won't propagate to the
+  // parameter, and only compute the cooling time.
+  OutputDustTemperature = -1;
+
+  Group_WriteAllData(DumpName, DumpNumber, TopGrid, MetaData, &Exterior
+#ifdef TRANSFER
+		       , ImplicitSolver
+#endif
+               );
+  
+
+  return SUCCESS;
+
+}

File src/enzo/PhotonGrid_Methods.h

View file
    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);
 			     float PhotonTestInitialFractionH2II,
 			     int RefineByOpticalDepth,
 			     int TotalRefinement,
-			     char *DensityFilename);
+			     char *DensityFilename,
+			     char *HIIFractionFilename,
+			     char *HeIIFractionFilename,
+			     char *HeIIIFractionFilename,
+			     char *TemperatureFilename);
 
 /************************************************************************/
 

File src/enzo/PhotonTestInitialize.C

View file
   char *dummy = new char[MAX_LINE_LENGTH];
   int   dim, ret, level, sphere, i, source;
   int   TotalRefinement;
+  dummy[0] = 0;
 
   /* set default parameters */
 
-  char *PhotonTestDensityFilename = NULL;
+  char *PhotonTestDensityFilename = NULL,
+    *PhotonTestHIIFractionFilename = NULL,
+    *PhotonTestHeIIFractionFilename = NULL,
+    *PhotonTestHeIIIFractionFilename = NULL,
+    *PhotonTestTemperatureFilename = NULL;
   int PhotonTestNumberOfSpheres = 1;
   int PhotonTestUseParticles    = FALSE;
   int PhotonTestUseColour       = FALSE;
 		  &PhotonTestInitialTemperature);
     if (sscanf(line, "PhotonTestDensityFilename = %s", dummy) == 1) {
       ret++;
-      PhotonTestDensityFilename = dummy;
+      PhotonTestDensityFilename = new char[MAX_LINE_LENGTH];
+      strcpy(PhotonTestDensityFilename, dummy);
+    }
+    if (sscanf(line, "PhotonTestHIIFractionFilename = %s", dummy) == 1) {
+      ret++;
+      PhotonTestHIIFractionFilename = new char[MAX_LINE_LENGTH];
+      strcpy(PhotonTestHIIFractionFilename, dummy);
+    }
+    if (sscanf(line, "PhotonTestHeIIFractionFilename = %s", dummy) == 1) {
+      ret++;
+      PhotonTestHeIIFractionFilename = new char[MAX_LINE_LENGTH];
+      strcpy(PhotonTestHeIIFractionFilename, dummy);
+    }
+    if (sscanf(line, "PhotonTestHeIIIFractionFilename = %s", dummy) == 1) {
+      ret++;
+      PhotonTestHeIIIFractionFilename = new char[MAX_LINE_LENGTH];
+      strcpy(PhotonTestHeIIIFractionFilename, dummy);
+    }
+    if (sscanf(line, "PhotonTestTemperatureFilename = %s", dummy) == 1) {
+      ret++;
+      PhotonTestTemperatureFilename = new char[MAX_LINE_LENGTH];
+      strcpy(PhotonTestTemperatureFilename, dummy);
     }
     ret += sscanf(line, "PhotonTestUniformVelocity = %"FSYM" %"FSYM" %"FSYM, 
 		  PhotonTestUniformVelocity, PhotonTestUniformVelocity+1,
 	     PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII,
 	     PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM,
 	     PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II, 
-	     RefineByOpticalDepth, TotalRefinement, PhotonTestDensityFilename);
+	     RefineByOpticalDepth, TotalRefinement, PhotonTestDensityFilename,
+	     PhotonTestHIIFractionFilename, PhotonTestHeIIFractionFilename,
+	     PhotonTestHeIIIFractionFilename, PhotonTestTemperatureFilename);
 
     CurrentGrid = CurrentGrid->NextGridThisLevel;
 
 	     PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII,
 	     PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM,
 	     PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II,
-	     RefineByOpticalDepth, TotalRefinement, NULL) == FAIL) {
+	     RefineByOpticalDepth, TotalRefinement, 
+	     NULL, NULL, NULL, NULL, NULL) == FAIL) {
 	  ENZO_FAIL("Error in PhotonTestInitializeGrid.\n");
 	}
 	Temp = Temp->NextGridThisLevel;
 		    PhotonTestInitialFractionHII, PhotonTestInitialFractionHeII,
 		    PhotonTestInitialFractionHeIII, PhotonTestInitialFractionHM,
 		    PhotonTestInitialFractionH2I, PhotonTestInitialFractionH2II,
-		    RefineByOpticalDepth, TotalRefinement, NULL) == FAIL) {
+		    RefineByOpticalDepth, TotalRefinement, 
+		    NULL, NULL, NULL, NULL, NULL) == FAIL) {
 	    ENZO_FAIL("Error in PhotonTestInitializeGrid.\n");
 	  }
 	  Temp = Temp->NextGridThisLevel;
   }
   } // ENDIF !Reinitialize
 
+  delete[] PhotonTestDensityFilename;
+  delete[] PhotonTestHIIFractionFilename;
+  delete[] PhotonTestHeIIFractionFilename;
+  delete[] PhotonTestHeIIIFractionFilename;
+  delete[] PhotonTestTemperatureFilename;
   delete [] dummy;
 
   return SUCCESS;

File src/enzo/ProjectToPlane2.C

View file
 	Units = "";
 	break;
       case 1:
-	Rho2Weighted = TRUE;
+	RhoWeighted = TRUE;
 	FieldName = "projected_rho2_weighted_temperature";
 	Units = "K";
 	break;
 	Units = "";
 	break;
       case 3:
-	Rho2Weighted = TRUE;
+	RhoWeighted = TRUE;
 	FieldName = "projected_metallicity";
 	Units = "relative to solar";
 	break;
 	Units = "";
 	break;
       case 5:
-	Rho2Weighted = TRUE;
+	RhoWeighted = TRUE;
 	FieldName = "projected_electron_fraction";
 	Units = "";
 	break;
       case 6:
-	Rho2Weighted = TRUE;
+	RhoWeighted = TRUE;
 	FieldName = "projected_H2_fraction";
 	Units = "";
 	break;

File src/enzo/RadiationFieldCalculateRates.C

View file
 
     aUnits = 1.0/(1.0 + InitialRedshift);
     Redshift = 1.0/(a*aUnits) - 1;
-  } else {  
+  } else if (RadiationFieldRedshift >= 0) {  
     Redshift = RadiationFieldRedshift;   
     CoolData.RadiationRedshiftOn = RadiationFieldRedshift+0.2;
     CoolData.RadiationRedshiftOff = 0.0;

File src/enzo/RadiationSource.h

File contents unchanged.

File src/enzo/RadiativeTransferParameters.h

View file
 
 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

View file
   RadiativeTransferPropagationDistance        = 0.1;
   RadiativeTransferCoupledRateSolver          = TRUE;
   RadiativeTransferOpticallyThinH2            = TRUE;
+  RadiativeTransferOpticallyThinXray          = FALSE;
   RadiativeTransferFluxBackgroundLimit        = 0.01;
   RadiativeTransferSplitPhotonRadius          = FLOAT_UNDEFINED; // kpc
   RadiativeTransferRaysPerCell                = 5.1;