Metal dust radiation

src/enzo/CallProblemSpecificRoutines.C

@@ -48,6 +48,7 @@ int CallProblemSpecificRoutines(TopGridData * MetaData, HierarchyEntry *ThisGrid
       FLOAT Center[] = {0,0,0}, CenterOfMass[3], DMCofM[3];
       ThisGrid->GridData->CalculateAngularMomentum
 	(Center, AM, MeanVelocity, DMVelocity, CenterOfMass, DMCofM);
+#ifdef TEMPORARY
       fprintf(stdout, 
 	      "level = %"ISYM" %"ISYM" %"ISYM"  "
 	      "Vel %"FSYM" %"FSYM" %"FSYM"  "
@@ -59,6 +60,7 @@ int CallProblemSpecificRoutines(TopGridData * MetaData, HierarchyEntry *ThisGrid
 	      DMVelocity[0], DMVelocity[1], DMVelocity[2],
 	      -CenterOfMass[0], -CenterOfMass[1], -CenterOfMass[2],
 	      DMCofM[0], DMCofM[1], DMCofM[2]);
+#endif
     }
     if (StellarWindSpeed > 0)
       ThisGrid->GridData->AddStellarWind();

src/enzo/CollapseTestInitialize.C

@@ -71,7 +71,44 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr,
   const char *DIName    = "DI_Density";
   const char *DIIName   = "DII_Density";
   const char *HDIName   = "HDI_Density";
+  const char *   HeHIIName =   "HeHII_Density";
+  const char *      DMName =      "DM_Density";
+  const char *    HDIIName =    "HDII_Density";
+  const char *      CIName =      "CI_Density";
+  const char *     CIIName =     "CII_Density";
+  const char *      COName =      "CO_Density";
+  const char *     CO2Name =     "CO2_Density";
+  const char *      OIName =      "OI_Density";
+  const char *      OHName =      "OH_Density";
+  const char *     H2OName =     "H2O_Density";
+  const char *      O2Name =      "O2_Density";
+  const char *     SiIName =     "SiI_Density";
+  const char *    SiOIName =    "SiOI_Density";
+  const char *   SiO2IName =   "SiO2I_Density";
+  const char *      CHName =      "CH_Density";
+  const char *     CH2Name =     "CH2_Density";
+  const char *    COIIName =    "COII_Density";
+  const char *     OIIName =     "OII_Density";
+  const char *    OHIIName =    "OHII_Density";
+  const char *   H2OIIName =   "H2OII_Density";
+  const char *   H3OIIName =   "H3OII_Density";
+  const char *    O2IIName =    "O2II_Density";
+  const char *      MgName =      "Mg_Density";
+  const char *      AlName =      "Al_Density";
+  const char *       SName =       "S_Density";
+  const char *      FeName =      "Fe_Density";
+  const char *     SiMName =     "SiM_Density";
+  const char *     FeMName =     "FeM_Density";
+  const char * Mg2SiO4Name = "Mg2SiO4_Density";
+  const char *  MgSiO3Name =  "MgSiO3_Density";
+  const char *   Fe3O4Name =   "Fe3O4_Density";
+  const char *      ACName =      "AC_Density";
+  const char *   SiO2DName =   "SiO2D_Density";
+  const char *     MgOName =     "MgO_Density";
+  const char *     FeSName =     "FeS_Density";
+  const char *   Al2O3Name =   "Al2O3_Density";
   const char *MetalName = "Metal_Density";
+  const char *DustName  = "Dust_Density";
 
   /* declarations */
 
@@ -557,11 +594,65 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr,
       DataLabel[count++] = (char*) DIIName;
       DataLabel[count++] = (char*) HDIName;
     }
+    if (MultiSpecies > 3) {
+      DataLabel[count++] = (char*) HeHIIName;
+      DataLabel[count++] = (char*) DMName;
+      DataLabel[count++] = (char*) HDIIName;
+    }
+    if (MetalChemistry > 0) {
+      DataLabel[count++] = (char*) CIName;
+      DataLabel[count++] = (char*) CIIName;
+      DataLabel[count++] = (char*) COName;
+      DataLabel[count++] = (char*) CO2Name;
+      DataLabel[count++] = (char*) OIName;
+      DataLabel[count++] = (char*) OHName;
+      DataLabel[count++] = (char*) H2OName;
+      DataLabel[count++] = (char*) O2Name;
+      DataLabel[count++] = (char*) SiIName;
+      DataLabel[count++] = (char*) SiOIName;
+      DataLabel[count++] = (char*) SiO2IName;
+      DataLabel[count++] = (char*) CHName;
+      DataLabel[count++] = (char*) CH2Name;
+      DataLabel[count++] = (char*) COIIName;
+      DataLabel[count++] = (char*) OIIName;
+      DataLabel[count++] = (char*) OHIIName;
+      DataLabel[count++] = (char*) H2OIIName;
+      DataLabel[count++] = (char*) H3OIIName;
+      DataLabel[count++] = (char*) O2IIName;
+      if (GrainGrowth || DustSublimation) {
+        if (DustSpecies > 0) {
+          DataLabel[count++] = (char*) MgName;
+        }
+        if (DustSpecies > 1) {
+          DataLabel[count++] = (char*) AlName;
+          DataLabel[count++] = (char*) SName;
+          DataLabel[count++] = (char*) FeName;
+        }
+      }
+    }
+    if (GrainGrowth || DustSublimation) {
+      if (DustSpecies > 0) {
+        DataLabel[count++] = (char*) MgSiO3Name;
+        DataLabel[count++] = (char*) ACName;
+      }
+      if (DustSpecies > 1) {
+        DataLabel[count++] = (char*) SiMName;
+        DataLabel[count++] = (char*) FeMName;
+        DataLabel[count++] = (char*) Mg2SiO4Name;
+        DataLabel[count++] = (char*) Fe3O4Name;
+        DataLabel[count++] = (char*) SiO2DName;
+        DataLabel[count++] = (char*) MgOName;
+        DataLabel[count++] = (char*) FeSName;
+        DataLabel[count++] = (char*) Al2O3Name;
+      }
+    }
   }  // if Multispecies
   if (CollapseTestUseColour)
     DataLabel[count++] = (char*) ColourName;
   if (CollapseTestUseMetals)
     DataLabel[count++] = (char*) MetalName;
+  if (UseDustDensityField)
+    DataLabel[count++] = (char*) DustName;
 
   for (i = 0; i < count; i++)
     DataUnits[i] = NULL;

src/enzo/CommunicationReceiverPhotons.C

@@ -198,6 +198,10 @@ int CommunicationReceiverPhotons(LevelHierarchyEntry *LevelArray[],
       NewPack->EmissionTime		= RecvBuffer[i].buffer.EmissionTime;
       NewPack->CurrentTime		= RecvBuffer[i].buffer.CurrentTime;
       NewPack->ColumnDensity	= RecvBuffer[i].buffer.ColumnDensity;
+      NewPack->ColumnDensityHDI = RecvBuffer[i].buffer.ColumnDensityHDI;
+      NewPack->ColumnDensityCO  = RecvBuffer[i].buffer.ColumnDensityCO;
+      NewPack->ColumnDensityOH  = RecvBuffer[i].buffer.ColumnDensityOH;
+      NewPack->ColumnDensityH2O = RecvBuffer[i].buffer.ColumnDensityH2O;
       NewPack->CrossSection		= RecvBuffer[i].buffer.CrossSection;
       NewPack->Radius		= RecvBuffer[i].buffer.Radius;
       NewPack->ipix			= RecvBuffer[i].buffer.ipix;

src/enzo/CommunicationTransferActiveParticles.C

@@ -92,7 +92,7 @@ int CommunicationTransferActiveParticles
 
     for (i = 0; i < Layout[dim]; i++) {
       ExactCount += ExactDims;
-      if (dim == 0)
+      if (i < Layout[dim]-1)
 	ThisCount = nint(0.5*ExactCount)*2 - DisplacementCount;
       else
 	ThisCount = nint(ExactCount) - DisplacementCount;

src/enzo/CommunicationTransferParticlesOpt.C

@@ -96,7 +96,7 @@ int CommunicationTransferParticles(grid *GridPointer[], int NumberOfGrids,
 
     for (i = 0; i < Layout[dim]; i++) {
       ExactCount += ExactDims;
-      if (dim == 0)
+      if (i < Layout[dim]-1)
 	ThisCount = nint(0.5*ExactCount)*2 - DisplacementCount;
       else
 	ThisCount = nint(ExactCount) - DisplacementCount;

src/enzo/CommunicationTransferPhotons.C

@@ -223,6 +223,14 @@ int CommunicationTransferPhotons(LevelHierarchyEntry *LevelArray[],
 	Mover->PhotonPackage->CurrentTime;
       SendList[ToProc][ToCount].buffer.ColumnDensity        = 
 	Mover->PhotonPackage->ColumnDensity;
+      SendList[ToProc][ToCount].buffer.ColumnDensityHDI     = 
+	Mover->PhotonPackage->ColumnDensityHDI;
+      SendList[ToProc][ToCount].buffer.ColumnDensityCO      = 
+	Mover->PhotonPackage->ColumnDensityCO;
+      SendList[ToProc][ToCount].buffer.ColumnDensityOH      = 
+	Mover->PhotonPackage->ColumnDensityOH;
+      SendList[ToProc][ToCount].buffer.ColumnDensityH2O     = 
+	Mover->PhotonPackage->ColumnDensityH2O;
       SendList[ToProc][ToCount].buffer.CrossSection         = 
 	Mover->PhotonPackage->CrossSection;
 

src/enzo/CommunicationTransferStarsOpt.C

@@ -94,7 +94,7 @@ int CommunicationTransferStars(grid *GridPointer[], int NumberOfGrids,
 
     for (i = 0; i < Layout[dim]; i++) {
       ExactCount += ExactDims;
-      if (dim == 0)
+      if (i < Layout[dim]-1)
 	ThisCount = nint(0.5*ExactCount)*2 - DisplacementCount;
       else
 	ThisCount = nint(ExactCount) - DisplacementCount;

src/enzo/CosmologySimulationInitialize.C

@@ -123,14 +123,53 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
   char *DIName    = "DI_Density";
   char *DIIName   = "DII_Density";
   char *HDIName   = "HDI_Density";
+  char *   HeHIIName =   "HeHII_Density";
+  char *      DMName =      "DM_Density";
+  char *    HDIIName =    "HDII_Density";
+  char *      CIName =      "CI_Density";
+  char *     CIIName =     "CII_Density";
+  char *      COName =      "CO_Density";
+  char *     CO2Name =     "CO2_Density";
+  char *      OIName =      "OI_Density";
+  char *      OHName =      "OH_Density";
+  char *     H2OName =     "H2O_Density";
+  char *      O2Name =      "O2_Density";
+  char *     SiIName =     "SiI_Density";
+  char *    SiOIName =    "SiOI_Density";
+  char *   SiO2IName =   "SiO2I_Density";
+  char *      CHName =      "CH_Density";
+  char *     CH2Name =     "CH2_Density";
+  char *    COIIName =    "COII_Density";
+  char *     OIIName =     "OII_Density";
+  char *    OHIIName =    "OHII_Density";
+  char *   H2OIIName =   "H2OII_Density";
+  char *   H3OIIName =   "H3OII_Density";
+  char *    O2IIName =    "O2II_Density";
+  char *      MgName =      "Mg_Density";
+  char *      AlName =      "Al_Density";
+  char *       SName =       "S_Density";
+  char *      FeName =      "Fe_Density";
+  char *     SiMName =     "SiM_Density";
+  char *     FeMName =     "FeM_Density";
+  char * Mg2SiO4Name = "Mg2SiO4_Density";
+  char *  MgSiO3Name =  "MgSiO3_Density";
+  char *   Fe3O4Name =   "Fe3O4_Density";
+  char *      ACName =      "AC_Density";
+  char *   SiO2DName =   "SiO2D_Density";
+  char *     MgOName =     "MgO_Density";
+  char *     FeSName =     "FeS_Density";
+  char *   Al2O3Name =   "Al2O3_Density";
   char *MetalName = "Metal_Density";
   char *MetalIaName = "MetalSNIa_Density";
+  char *DustName = "Dust_Density";
   char *GPotName  = "Grav_Potential";
   char *ForbidName  = "ForbiddenRefinement";
   char *MachName   = "Mach";
   char *PSTempName = "PreShock_Temperature";
   char *PSDenName  = "PreShock_Density";
-  char *ExtraNames[2] = {"Z_Field1", "Z_Field2"};
+  char *ExtraNames[12] = {"Z_Field0", "Z_Field1", "Z_Field2", "Z_Field3"
+                        , "Z_Field4", "Z_Field5", "Z_Field6", "Z_Field7"
+                        , "Z_Field8", "Z_Field9", "Z_Field10","Z_Field11"};
   char *BxName = "Bx";
   char *ByName = "By";
   char *BzName = "Bz";
@@ -786,6 +825,58 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
       DataLabel[i++] = DIIName;
       DataLabel[i++] = HDIName;
     }
+    if (MultiSpecies > 3) {
+      DataLabel[i++] = HeHIIName;
+      DataLabel[i++] = DMName;
+      DataLabel[i++] = HDIIName;
+    }
+    if (MetalChemistry > 0) {
+      DataLabel[i++] = CIName;
+      DataLabel[i++] = CIIName;
+      DataLabel[i++] = COName;
+      DataLabel[i++] = CO2Name;
+      DataLabel[i++] = OIName;
+      DataLabel[i++] = OHName;
+      DataLabel[i++] = H2OName;
+      DataLabel[i++] = O2Name;
+      DataLabel[i++] = SiIName;
+      DataLabel[i++] = SiOIName;
+      DataLabel[i++] = SiO2IName;
+      DataLabel[i++] = CHName;
+      DataLabel[i++] = CH2Name;
+      DataLabel[i++] = COIIName;
+      DataLabel[i++] = OIIName;
+      DataLabel[i++] = OHIIName;
+      DataLabel[i++] = H2OIIName;
+      DataLabel[i++] = H3OIIName;
+      DataLabel[i++] = O2IIName;
+      if (GrainGrowth || DustSublimation) {
+        if (DustSpecies > 0) {
+          DataLabel[i++] = MgName;
+        }
+        if (DustSpecies > 1) {
+          DataLabel[i++] = AlName;
+          DataLabel[i++] = SName;
+          DataLabel[i++] = FeName;
+        }
+      }
+    }
+    if (GrainGrowth || DustSublimation) {
+      if (DustSpecies > 0) {
+        DataLabel[i++] = MgSiO3Name;
+        DataLabel[i++] = ACName;
+      }
+      if (DustSpecies > 1) {
+        DataLabel[i++] = SiMName;
+        DataLabel[i++] = FeMName;
+        DataLabel[i++] = Mg2SiO4Name;
+        DataLabel[i++] = Fe3O4Name;
+        DataLabel[i++] = SiO2DName;
+        DataLabel[i++] = MgOName;
+        DataLabel[i++] = FeSName;
+        DataLabel[i++] = Al2O3Name;
+      }
+    }
   }
   if (CosmologySimulationUseMetallicityField) {
     DataLabel[i++] = MetalName;
@@ -794,8 +885,23 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
     if(MultiMetals){
       DataLabel[i++] = ExtraNames[0];
       DataLabel[i++] = ExtraNames[1];
+      DataLabel[i++] = ExtraNames[2];
+      DataLabel[i++] = ExtraNames[3];
+      DataLabel[i++] = ExtraNames[4];
+      DataLabel[i++] = ExtraNames[5];
+      DataLabel[i++] = ExtraNames[6];
+      DataLabel[i++] = ExtraNames[7];
+      DataLabel[i++] = ExtraNames[8];
+      DataLabel[i++] = ExtraNames[9];
+      DataLabel[i++] = ExtraNames[10];
+      DataLabel[i++] = ExtraNames[11];
     }
   }
+
+  if (UseDustDensityField) {
+    DataLabel[i++] = DustName;
+  }
+
   if(STARMAKE_METHOD(COLORED_POP3_STAR)){
     DataLabel[i++] = ForbidName;
   }

src/enzo/EvolvePhotons.C

@@ -303,7 +303,7 @@ int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 	    ENZO_FAIL("Error in InitializeTemperatureFieldForComptonHeating.\n");
 	  }	
     /* Initialize Temperature Field for H2 shielding approximation */
-    if(RadiativeTransferH2ShieldType == 1 || ProblemType == 50) {
+    if(RadiativeTransferH2ShieldType == 1 || (ProblemType == 50 && RadiativeTransferUseH2Shielding)) {
       for (lvl = MAX_DEPTH_OF_HIERARCHY-1; lvl >= 0 ; lvl--)
 	for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel) 
 	  if (Temp->GridData->InitializeTemperatureFieldForH2Shield() == FAIL) {
@@ -639,7 +639,7 @@ int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 	  if (Temp->GridData->FinalizeTemperatureFieldForComptonHeating() == FAIL) {  
 	    ENZO_FAIL("Error in FinalizeTemperatureFieldForComptonHeating.\n");
 	  }	
-    if (RadiativeTransferH2ShieldType == 1 || ProblemType == 50)
+    if(RadiativeTransferH2ShieldType == 1 || (ProblemType == 50 && RadiativeTransferUseH2Shielding))
       for (lvl = 0; lvl < MAX_DEPTH_OF_HIERARCHY-1; lvl++)
 	for (Temp = LevelArray[lvl]; Temp; Temp = Temp->NextGridThisLevel)
 	  if (Temp->GridData->FinalizeTemperatureFieldForH2Shield() == FAIL) {