Import of Shock/Cosmic Ray Code. Modifies quite a few files, but does…

…n't affect anything unless CRModel > 0.

--HG--
branch : week-of-code

src/enzo/CosmologySimulationInitialize.C

@@ -111,6 +111,10 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
   char *HDIName   = "HDI_Density";
   char *MetalName = "Metal_Density";
   char *GPotName  = "Grav_Potential";
+  char *MachName   = "Mach";
+  char *CRName     = "CR_Density";
+  char *PSTempName = "PreShock_Temperature";
+  char *PSDenName  = "PreShock_Density";
   char *ExtraNames[2] = {"Z_Field1", "Z_Field2"};
 
 
@@ -589,7 +593,14 @@ int CosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
       DataLabel[i++] = HDIName;
     }
   }
-
+  if (CRModel) {
+    DataLabel[i++] = MachName;
+    if(StorePreShockFields){
+      DataLabel[i++] = PSTempName;
+      DataLabel[i++] = PSDenName;
+    }
+    DataLabel[i++] = CRName;
+  } 
   if (CosmologySimulationUseMetallicityField) {
     DataLabel[i++] = MetalName;
     if(MultiMetals){

src/enzo/EvolveLevel.C

@@ -34,6 +34,8 @@
 /                Optional StaticSiblingList for root grid
 /  modified8:  April, 2009 by John Wise
 /                Added star particle class and radiative transfer
+/  modified9:  July, 2009 by Sam Skillman
+/                Added shock and cosmic ray analysis
 /
 /  PURPOSE:
 /    This routine is the main grid evolution function.  It assumes that the
@@ -417,6 +419,10 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[],
 
       Grids[grid1]->GridData->StarParticleHandler
 	(Grids[grid1]->NextGridNextLevel, level);
+
+      /* Include shock-finding and cosmic ray acceleration */
+
+      Grids[grid1]->GridData->ShocksHandler();
 
       /* Gravity: clean up AccelerationField. */
 

src/enzo/Grid.h

@@ -512,6 +512,10 @@ class grid
 
    int MultiSpeciesHandler();
 
+/* Handle the selection of shock finding algorithm */
+
+   int ShocksHandler();
+
 /* Solve the radiative cooling/heating equations  */
 
    int SolveRadiativeCooling();
@@ -528,6 +532,14 @@ class grid
 
    int RadiationComputeDensities(int level);
 
+// -------------------------------------------------------------------------
+// Functions for shock finding and cosmic ray acceleration
+//
+   int FindShocks();
+   int FindTempSplitShocks();
+   int FindVelShocks();
+   int FindVelSplitShocks();
+
 // -------------------------------------------------------------------------
 // Functions for grid (re)generation.
 //
@@ -1220,6 +1232,10 @@ void SortParticlesByNumber();
 			    int &HMNum, int &H2INum, int &H2IINum,
                             int &DINum, int &DIINum, int &HDINum);
 
+/* Identify shock/cosmic ray fields. */
+  int IdentifyCRSpeciesFields(int &MachNum, int&CRNum, 
+			      int &PSTempNum, int &PSDenNum);
+
 /* Zeus Solver. */
 
   int ZeusSolver(float *gamma, int igamfield, int nhy, 
@@ -1262,6 +1278,9 @@ void SortParticlesByNumber();
 // Functions for Specific problems (usually problem generator functions).
 //
 
+/* Generalized Extra Field Grid Initializer (returns NumberOfBaryonFields) */
+  int InitializeTestProblemGrid(int field_counter);
+
 /* Protostellar Collapse problem: initialize grid (returns SUCCESS or FAIL) */
   int ProtostellarCollapseInitializeGrid(float CoreDensity,
 					 float CoreEnergy,

src/enzo/Grid_CorrectForRefinedFluxes.C

@@ -235,8 +235,9 @@ int grid::CorrectForRefinedFluxes(fluxes *InitialFluxes,
 	    total number density summed over ionization, etc.) */
 
 	for (field = 0; field < NumberOfBaryonFields; field++)
-	  if ((RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
-		 			 FieldType[field] != InternalEnergy))
+ 	  if (FieldTypeNoInterpolate(FieldType[field]) == FALSE &&
+ 	      (RadiativeCooling == 0 || (FieldType[field] != TotalEnergy &&
+					 FieldType[field] != InternalEnergy))
 	      && (FieldType[field] < ElectronDensity)) {
 	  for (k = Start[2]; k <= End[2]; k++)
 	    for (j = Start[1]; j <= End[1]; j++)

src/enzo/Grid_CosmologySimulationInitializeGrid.C

@@ -117,6 +117,8 @@ int grid::CosmologySimulationInitializeGrid(
   int idim, dim, i, j, vel, OneComponentPerFile, ndim, level;
   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
       DINum, DIINum, HDINum, MetalNum;
+
+  int CRNum, MachNum, PSTempNum, PSDenNum;
 
   int ExtraField[2];
 
@@ -250,6 +252,21 @@ int grid::CosmologySimulationInitializeGrid(
     }
     if (WritePotential)
       FieldType[NumberOfBaryonFields++] = GravPotential;
+    if(CRModel){
+      FieldType[MachNum   = NumberOfBaryonFields++] = Mach;
+      if(StorePreShockFields){
+	FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature;
+	FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity;
+      }
+      FieldType[CRNum     = NumberOfBaryonFields++] = CRDensity;
+    }
+    if (UseMetallicityField) {
+      FieldType[MetalNum = NumberOfBaryonFields++] = Metallicity;
+      if(MultiMetals){
+	FieldType[ExtraField[0] = NumberOfBaryonFields++] = ExtraType0;
+	FieldType[ExtraField[1] = NumberOfBaryonFields++] = ExtraType1;
+      }
+    }
   }
 
   // Set the subgrid static flag
@@ -421,8 +438,20 @@ int grid::CosmologySimulationInitializeGrid(
 	                             BaryonField[H2INum][i];
       }
 
+      //Shock/Cosmic Ray Model
+      if(CRModel && ReadData)
+	for (i = 0; i < size; i++){
+	  BaryonField[MachNum][i] = tiny_number;
+	  BaryonField[CRNum][i] = tiny_number;
+	  if(StorePreShockFields){
+	    BaryonField[PSTempNum][i] = tiny_number;
+	    BaryonField[PSDenNum][i] = tiny_number;
+	  }
+
+	}
+
     }
-
+  
   // If using metallicity, set the field
 
   if (UseMetallicityField && ReadData)

src/enzo/Grid_InterpolateBoundaryFromParent.C

@@ -361,7 +361,7 @@ int grid::InterpolateBoundaryFromParent(grid *ParentGrid)
 
       if (FieldType[field] == Density)
 	FieldPointer = TemporaryDensityField;
-      else
+      else if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
 	FieldPointer = TemporaryField;
 
       /* Copy needed portion of temp field to current grid. */

src/enzo/Grid_InterpolateFieldValues.C

@@ -327,8 +327,8 @@ int grid::InterpolateFieldValues(grid *ParentGrid)
 
       if (FieldType[field] == Density)
 	FieldPointer = TemporaryDensityField;
-      else
-	FieldPointer = TemporaryField;
+      else if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
+	  FieldPointer = TemporaryField;
 
       /* Copy needed portion of temp field to current grid. */
 

src/enzo/Grid_ProjectSolutionToParentGrid.C

@@ -120,9 +120,10 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
 
   if (ProcessorNumber == MyProcessorNumber)
     for (field = 0; field < NumberOfBaryonFields; field++)
-      if (FieldTypeIsDensity(FieldType[field]) == FALSE && (
-	  (FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
-          || HydroMethod != Zeus_Hydro      )       )
+      if ((FieldTypeIsDensity(FieldType[field]) == FALSE &&  
+	   ((FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
+	    || HydroMethod != Zeus_Hydro      )       ) &&
+	  (FieldTypeNoInterpolate(FieldType[field]) == FALSE))
       FORTRAN_NAME(mult3d)(BaryonField[DensNum], BaryonField[field],
 			   &Size, &One, &One, &Size, &One, &One,
 			   &Zero, &Zero, &Zero, &Zero, &Zero, &Zero);
@@ -184,8 +185,9 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
 
 	  for (i = 0; i < Dim[0]; i += skipi) {
 	    i1 = i/Refinement[0];
-	    ParentGrid.BaryonField[field][pindex+i1] +=
-	      BaryonField[field][gindex+i]*weight;
+	    if (FieldTypeNoInterpolate(FieldType[field]) == FALSE)
+	      ParentGrid.BaryonField[field][pindex+i1] +=
+		BaryonField[field][gindex+i]*weight;
 	  }
 	}
       }
@@ -222,9 +224,10 @@ int grid::ProjectSolutionToParentGrid(grid &ParentGrid)
   /* Divide all fields by mass to return to original quantity. */
 
   for (field = 0; field < NumberOfBaryonFields; field++)
-    if (FieldTypeIsDensity(FieldType[field]) == FALSE && (
+    if ((FieldTypeIsDensity(FieldType[field]) == FALSE && (
 	  (FieldType[field] < Velocity1 || FieldType[field] > Velocity3)
-          || HydroMethod != Zeus_Hydro      )       ) {
+          || HydroMethod != Zeus_Hydro      )       ) &&
+	(FieldTypeNoInterpolate(FieldType[field]) == FALSE)){
       if (ProcessorNumber == MyProcessorNumber)
 	FORTRAN_NAME(div3d)(BaryonField[DensNum], BaryonField[field],
 			    &Size, &One, &One, &Size, &One, &One,

src/enzo/Grid_SolveHydroEquations.C

@@ -121,8 +121,44 @@ int grid::SolveHydroEquations(int CycleNumber, int NumberOfSubgrids,
       colnum[NumberOfColours++] = SNColourNum;
     }
 
-    /* Determine if Gamma should be a scalar or a field. */
+    /* Add shock/cosmic ray variables as a colour variable. */
 
+    if(CRModel){
+      int MachNum, CRNum, PSTempNum,PSDenNum;
+
+      if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) {
+	ENZO_FAIL("Error in IdentifyCRSpeciesFields.")
+      }
+
+      colnum[NumberOfColours++] = MachNum;
+      if(StorePreShockFields){
+	colnum[NumberOfColours++] = PSTempNum;
+	colnum[NumberOfColours++] = PSDenNum;
+      }
+      colnum[NumberOfColours++] = CRNum;
+
+      /*  Zero out Mach number array so that we don't get strange advection */
+
+      float *mach = BaryonField[MachNum];
+      if(StorePreShockFields){
+	float *pstemp = BaryonField[PSTempNum];
+	float *psden = BaryonField[PSDenNum];
+	for (i = 0; i < size; i++){
+	  pstemp[i] = tiny_number;
+	  psden[i] = tiny_number;
+	}
+      }
+      for (i = 0; i < size; i++)
+	mach[i] = tiny_number;
+      if(CRModel == 2){  //zero out CR to get instantanous injection
+	float *cr = BaryonField[CRNum];
+	for (i = 0; i < size; i++)
+	  cr[i] = tiny_number;
+      }
+    }
+
+    /* Determine if Gamma should be a scalar or a field. */
+
     int UseGammaField = FALSE;
     float *GammaField;
     if (HydroMethod == Zeus_Hydro && MultiSpecies > 1) {

src/enzo/Make.config.objects

@@ -258,6 +258,7 @@ OBJS_CONFIG_LIB = \
         Grid_destructor.o \
         Grid_DetachForcingFromBaryonFields.o \
         Grid_DoubleMachInitializeGrid.o \
+	Grid_FindShocks.o \
         Grid_FastSiblingLocatorAddGrid.o \
         Grid_FastSiblingLocatorFindSiblings.o \
         Grid_FindAllStarParticles.o \
@@ -282,6 +283,7 @@ OBJS_CONFIG_LIB = \
         Grid_Group_ReadGrid.o \
         Grid_Group_WriteGridInterpolate.o \
         Grid_Group_WriteGrid.o \
+	Grid_IdentifyCRSpeciesFields.o \
         Grid_IdentifyNewSubgrids.o \
         Grid_IdentifyNewSubgridsSmall.o \
         Grid_IdentifyPhysicalQuantities.o \
@@ -350,6 +352,7 @@ OBJS_CONFIG_LIB = \
         Grid_SetParticleMassFlaggingField.o \
         Grid_ShearingBoxInitializeGrid.o \
         Grid_ShockTubeInitializeGrid.o \
+	Grid_ShocksHandler.o \
         Grid_SolveForPotential.o \
         Grid_SolveHydroEquations.o \
         Grid_SolveRadiativeCooling.o \
@@ -369,6 +372,7 @@ OBJS_CONFIG_LIB = \
         Grid_TracerParticleCreateParticles.o \
         Grid_TracerParticleOutputData.o \
         Grid_TracerParticleSetVelocity.o \
+	Grid_TestProblemInitializeGrid.o \
         Grid_TransferSubgridParticles.o \
         Grid_TransferSubgridStars.o \
         Grid_TurbulenceSimulationInitialize.o \
@@ -399,6 +403,7 @@ OBJS_CONFIG_LIB = \
         IdentifyNewSubgridsBySignature.o \
         ImplosionInitialize.o \
         InitializeCloudyCooling.o \
+	InitializeCosmicRayData.o \
         InitializeEquilibriumCoolData.o \
         InitializeLocal.o \
         InitializeMovieFile.o \
@@ -558,6 +563,7 @@ OBJS_CONFIG_LIB = \
         TestGravitySphereInitialize.o \
         TestOrbitInitialize.o \
         TracerParticleCreation.o \
+	TestProblemDataInitialize.o \		
         tscint1d.o \
         tscint2d.o \
         tscint3d.o \