modified star_maker7(INSTANT_STAR) so that it wouldn't create stars w…

…hen a MBH particle exists in that cell

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

input/make_Zcool_table.pro

@@ -9,7 +9,7 @@ pop3_on = 0
 
 h2_fraction = 1e-5
 density = 1.0
-redshift = 15.0
+redshift = 0.0 ;15.0
 h_fraction = 0.76  ;; Hydrogen mass fraction
 metallicity = 0.1  ;; relative to solar (don't matter here ... everything 
                    ;; approx. scales with Z)
@@ -39,7 +39,7 @@ printf, 1, "# ", t_bins, ne_bins
 printf, 1, "# ", t_range[0], t_range[1], ne_range[0], ne_range[1]
 printf, 1, "# Rates are in erg/s * cm^3 for [Z/H] = ", $
         STRING(ALOG10(metallicity), FORMAT='(F0.2)')
-printf, 1, "# Rates assumed to scale linearly with metallicity in the code."
+printf, 1, "# Rates assumed to scale linearly with metallicity in enzo (see cool1d_multi.src)."
 printf, 1, "#"
 line_format = STRING(ne_bins, FORMAT='("(",I0,"(G12.5))")')
 printf, 1, TRANSPOSE(total_rate), FORMAT=line_format

src/enzo/Grid_ParticleSplitter.C

@@ -32,7 +32,7 @@
 #include "StarParticleData.h"
 
 #define NO_DEBUG_PS 
-#define PARTICLE_IN_GRID_CHECK 
+#define PARTICLE_IN_GRID_CHECK //#####
 
 /* function prototypes */
 

src/enzo/Grid_StarParticleHandler.C

@@ -31,7 +31,7 @@
 
 #define  PROTONMASS  1.6726e-24
 
-#define PARTICLE_IN_GRID_CHECK  
+#define PARTICLE_IN_GRID_CHECK   //#####
 
 /* function prototypes */
 
@@ -110,9 +110,9 @@ extern "C" void FORTRAN_NAME(star_maker4)(int *nx, int *ny, int *nz,
      		 int *ibuff,
              int *imetal, hydro_method *imethod, float *mintdyn,
              float *odthresh, float *massff, float *smthrest, int *level,
-		 int *np,
+		 int *np, int *npart,
              FLOAT *xp, FLOAT *yp, FLOAT *zp, float *up, float *vp, float *wp,
-             float *mp, float *tdp, float *tcp, float *metalf);
+             float *mp, float *tdp, float *tcp, float *metalf, int *type, int *ctype);
 
 extern "C" void FORTRAN_NAME(star_maker5)(int *nx, int *ny, int *nz,
              float *d, float *dm, float *temp, float *coolrate, float *u, float *v, float *w,
@@ -379,9 +379,10 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
   /* Set variables to type defines to pass to FORTRAN routines */
 
   int NormalStarType = PARTICLE_TYPE_STAR;
+  int SinkParticleType = PARTICLE_TYPE_MUST_REFINE;
   int SingleStarType = PARTICLE_TYPE_SINGLE_STAR;
   int StarClusterType = PARTICLE_TYPE_CLUSTER;
-  int SinkParticleType = PARTICLE_TYPE_MUST_REFINE;
+  int MBHParticleType = PARTICLE_TYPE_MBH;
   int ColorStar = PARTICLE_TYPE_COLOR_STAR;
 
   /* Compute the redshift. */
@@ -625,6 +626,11 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
     }
 
     if (STARMAKE_METHOD(COLORED_POP3_STAR)) {
+
+      //---- POPULATION III (SINGLE STAR)
+
+      NumberOfNewParticlesSoFar = NumberOfNewParticles;
+
       PFORTRAN_NAME(pop3_color_maker)
         (GridDimension, GridDimension+1, GridDimension+2, 
          BaryonField[DensNum], dmfield, 
@@ -637,7 +643,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
          &HydroMethod, &PopIIIColorDensityThreshold, 
            &level, &NumberOfNewParticles, 
          tg->ParticlePosition[0], tg->ParticlePosition[1], tg->ParticlePosition[2],
-          tg->ParticleVelocity[0], tg->ParticleVelocity[1], tg->ParticleVelocity[2],
+         tg->ParticleVelocity[0], tg->ParticleVelocity[1], tg->ParticleVelocity[2],
          tg->ParticleMass, tg->ParticleAttribute[2], tg->ParticleType, &ColorStar);
 
     }
@@ -692,7 +698,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
 	  JeansLengthRefinement = RefineByJeansLengthSafetyFactor;
       }
 
-      int MBHParticleType = -PARTICLE_TYPE_MBH; //minus sign needed for yet-to-be born Star particle
+      MBHParticleType = -PARTICLE_TYPE_MBH; //minus sign needed for yet-to-be born Star particle
 
       NumberOfNewParticlesSoFar = NumberOfNewParticles;
 
@@ -724,7 +730,7 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
     if (STARMAKE_METHOD(INSTANT_STAR) && level == MaximumRefinementLevel) {
 
       //---- MODIFIED SF ALGORITHM (NO-JEANS MASS, NO dt DEPENDENCE, NO STOCHASTIC SF, 
-      //                            only at MaximumRefinementLevel)
+      //                            only at MaximumRefinementLevel, avoids SF when MBH exists)
 
       NumberOfNewParticlesSoFar = NumberOfNewParticles;
 
@@ -739,13 +745,13 @@ int grid::StarParticleHandler(HierarchyEntry* SubgridPointer, int level)
           CellLeftEdge[2], &GhostZones,
        &MetallicityField, &HydroMethod, &StarMakerMinimumDynamicalTime,
        &StarMakerOverDensityThreshold, &StarMakerMassEfficiency,
-       &StarMakerMinimumMass, &level, &NumberOfNewParticles, 
+       &StarMakerMinimumMass, &level, &NumberOfNewParticles, &NumberOfParticles,
        tg->ParticlePosition[0], tg->ParticlePosition[1],
           tg->ParticlePosition[2],
        tg->ParticleVelocity[0], tg->ParticleVelocity[1],
           tg->ParticleVelocity[2],
        tg->ParticleMass, tg->ParticleAttribute[1], tg->ParticleAttribute[0],
-          tg->ParticleAttribute[2]);
+          tg->ParticleAttribute[2], ParticleType, &MBHParticleType);
 
       for (i = NumberOfNewParticlesSoFar; i < NumberOfNewParticles; i++)
           tg->ParticleType[i] = NormalStarType;

src/enzo/Make.config.objects

@@ -129,6 +129,7 @@ OBJS_CONFIG_LIB = \
         ExternalBoundary_AppendForcingToBaryonFields.o \
         ExternalBoundary_constructor.o \
         ExternalBoundary_DetachForcingFromBaryonFields.o \
+        ExternalBoundary_DeleteObsoleteFields.o \
         ExternalBoundary_IdentifyPhysicalQuantities.o \
         ExternalBoundary_InitializeExternalBoundaryFaceIO.o \
         ExternalBoundary_Prepare.o \
@@ -277,6 +278,7 @@ OBJS_CONFIG_LIB = \
         Grid_DeleteAllButParticles.o \
         Grid_DeleteAllFields.o \
         Grid_DeleteFlaggingField.o \
+        Grid_DeleteObsoleteFields.o \
 	Grid_FlagCellsToBeRefinedByMetallicity.o	\
         Grid_DepositMustRefineParticles.o \
 	Grid_FlagCellsToBeRefinedByMustRefineRegion.o	\

src/enzo/Makefile

@@ -47,6 +47,7 @@ MODULES  =
 VERBOSE  = 0
 
 # please remove this
+#SVN      = hg
 
 #-----------------------------------------------------------------------
 # Make.config.settings is used for setting default values to all compile-time 

src/enzo/ReadMetalCoolingRates.C

@@ -81,7 +81,6 @@ int ReadMetalCoolingRates(float TemperatureUnits, float LengthUnits,
 
   int prev_pos;
   int ixe, itemp, index, nbins;
-
   nbins = CoolData.NumberOfTemperatureBins * CoolData.NumberOfElectronFracBins;
   CoolData.metals = new float[nbins];
 
@@ -122,6 +121,28 @@ int ReadMetalCoolingRates(float TemperatureUnits, float LengthUnits,
   for (i = 0; i < nbins; i++)
     CoolData.metals[i] /= coolunit;
 
+#define OUTPUT_METAL_COOLING
+
+#ifdef OUTPUT_METAL_COOLING
+  /* Output cooling rate in code units. */
+
+  FILE *fptr2;
+
+  fptr2 = fopen("metal_cool_rates.out", "w");
+  for (itemp = 0; itemp < CoolData.NumberOfTemperatureBins; itemp++) {
+    fprintf(fptr2, "%"GSYM"\t", log10(CoolData.TemperatureStart) +
+	   (log10(CoolData.TemperatureEnd)-
+	    log10(CoolData.TemperatureStart)) * float(itemp)/
+	    float(CoolData.NumberOfTemperatureBins-1));
+    for (ixe = 0; ixe < CoolData.NumberOfElectronFracBins; ixe++) {
+      index = ixe*CoolData.NumberOfTemperatureBins + itemp;
+      fprintf(fptr2, "%"GSYM"\t", CoolData.metals[index]);
+    }
+    fprintf(fptr2, "\n");
+  }
+  fclose(fptr2);
+#endif
+
   return SUCCESS;
 
 }

src/enzo/WriteParameterFile.C

@@ -353,6 +353,7 @@ int WriteParameterFile(FILE *fptr, TopGridData &MetaData)
   fprintf(fptr, "RadiationFieldLevelRecompute   = %"ISYM"\n", RadiationFieldLevelRecompute);
   fprintf(fptr, "RadiationSpectrumNormalization = %"GSYM"\n", CoolData.f3);
   fprintf(fptr, "RadiationSpectrumSlope         = %"GSYM"\n", CoolData.alpha0);
+  fprintf(fptr, "PhotoelectricHeating           = %lf\n", PhotoelectricHeating*uheat);
 
   if (CoolData.ParameterFilename != NULL)
     fprintf(fptr, "CoolDataParameterFile = %s\n\n", CoolData.ParameterFilename);

src/enzo/cool1d_multi.src

@@ -684,13 +684,34 @@ c
      $          metala(indixe(i),   xe_idx+1) * (1-dlogxe) * (dlogT)   +
      $          metala(indixe(i)+1, xe_idx+1) * (dlogxe  ) * (dlogT)
 c
+c---------------------------------------------------------------------
+c #####
+c
+c               write(*,*) "------------"
+c               write(*,*) "tgas(i) = ", tgas(i)
+c               write(*,*) "tgasold(i) = ", tgasold(i)
+c               write(*,*) "logtem(i) = ", logtem(i)
+c               write(*,*) "dlogT = ", dlogT
+c               write(*,*) "indixe(i) = ", indixe(i)
+c               write(*,*) "--"
+c               write(*,*) "xe = ", xe
+c               write(*,*) "log_xe = ", log_xe
+c               write(*,*) "dlogxe = ", dlogxe
+c               write(*,*) "xe_idx = ", xe_idx
+c               write(*,*) "--"
+c               write(*,*) "metala(indixe(i), xe_idx) = ", 
+c     $                     metala(indixe(i), xe_idx)
+c
+c
+c---------------------------------------------------------------------
+c
 c  Subtract metal cooling at T_cmb
 c
             metalc_cmb = 
      $        metala(tcmb_idx,   xe_idx  ) * (1-dlogxe) * (1-dlogTcmb) +
      $        metala(tcmb_idx+1, xe_idx  ) * (dlogxe  ) * (1-dlogTcmb) +
      $        metala(tcmb_idx,   xe_idx+1) * (1-dlogxe) * (dlogTcmb)   +
-     $          metala(tcmb_idx+1, xe_idx+1) * (dlogxe  ) * (dlogTcmb)
+     $        metala(tcmb_idx+1, xe_idx+1) * (dlogxe  ) * (dlogTcmb)
                metalc(i) = metalc(i) - metalc_cmb
 
                xi = min(3.0, metal(i,j,k)/(d(i,j,k)*ZSOLAR))

src/enzo/star_maker7.src

@@ -10,9 +10,10 @@ c
      &                      d1, x1, v1, t1,
      &                      nmax, xstart, ystart, zstart, ibuff, 
      &                      imetal, imethod, mintdyn,
-     &                      odthresh, masseff, smthresh, level, np,
+     &                      odthresh, masseff, smthresh, level, 
+     &                      np, npart,
      &                      xp, yp, zp, up, vp, wp,
-     &                      mp, tdp, tcp, metalf)
+     &                      mp, tdp, tcp, metalf, type, ctype)
 c
 c  CREATES GALAXY PARTICLES
 c
@@ -34,8 +35,9 @@ C    imported star_maker2.src as of 12 June 2009 and changed the followings:
 c    [1] removed Jeans mass criteria (just like star_maker3.src)
 c    [2] removed stochastic star formation (unlike star_maker3.src)
 c    [3] removed dt dependence in stellar mass created 
-c    
-cc
+c    [4] prevents star formation when there is a MBH particle in the cell
+c        MBH particle will suck out the gas later in StarParticleFinalize
+c
 c
 c  INPUTS:
 c
@@ -63,7 +65,8 @@ c    procnum - processor number (for output)
 c
 c  OUTPUTS:
 c
-c    np   - number of particles created
+c    np   - number of particles created (=NumberOfNewParticles)
+c    npart   - number of particles already in here (=NumberOfParticles)
 c    x/y/z start - starting position of grid origin
 c    xp,yp,zp - positions of created particles
 c    up,vp,wp - velocities of created particles
@@ -80,8 +83,9 @@ c-----------------------------------------------------------------------
 c
 c  Arguments
 c
-      integer nx, ny, nz, ibuff, nmax, np, level, imetal, imethod
-      integer procnum
+      integer nx, ny, nz, ibuff, nmax, np, npart, level, imetal, imethod
+      integer procnum, ctype
+      integer type(npart)
       real    d(nx,ny,nz), dm(nx,ny,nz), temp(nx,ny,nz)
       real    u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz)
       real    r(nx,ny,nz), cooltime(nx,ny,nz), metal(nx,ny,nz)
@@ -99,7 +103,7 @@ c
 c
 c  Locals:
 c
-      integer  i, j, k, ii
+      integer  i, j, k, ii, n
       real   div, tdyn, dtot
       real   pi, G, sndspdC
       real   isosndsp2, starmass, starfraction, bmass, jeanmass
@@ -109,6 +113,16 @@ c
 c
       ii = np
 c
+c  check whether we have MBH particle here, in which case we skip the
+c  star formation routine
+c
+      do n=1, npart
+         if (type(n) .eq. ctype) then
+c            write(6,*) 'star_maker7: the cell has a MBH particle; move on'
+            goto 20
+         endif
+      enddo
+c
 c  for each zone, : "star" particle is created if answers to all the
 c  following questions are affirmative:
 c
@@ -263,7 +277,7 @@ c
 c      if (np .ne. 0) then
 c         write(6,*) 'Stars created: number,time,level: ', np, t, level
 c      endif
-cc
+c
       return
       end
 c