Stephen Skory’s Research

because it turns out using this multi-user wordpress kind of sucks.

I have done this to my enzohop today:

• Added a function to read in the cosmology parameters from the restart file so masses get calculated correctly.
• Added a command line parameter for HOP density threshold.
• Changed the default cosmology parameters.

• I added a page describing how to use the merger tree code on the yt site. I added functionality to choose which haloes at the end to inspect, and I fixed my box problems which had to do with the range I was instantiating the snapshots. Matt has committed the code to the trunk and is making some changes to how it fits into yt.

I’ve put the Kraken work on SDSC HPSS as skory_kraken.tar

This is my best guess at what’s going, which leads to a NaN in star_maker5.src. Take a look at the diagram below, which isn’t a complete picture of what’s going on, but was enough for me to decide I dislike Gadget.

The way to read the graph is to realize I started at the top, in star_maker5.src. I saw that pstar was NaN. From that point on I’m mostly tracking zeros. The problems are identified along each arrow. So, pstar is NaN because x=NaN. x=NaN because in the expression shown, y=0. And so on. Each grey box represents a different source file.

The final boxes are on the right side, where I tried to hack this stuff by setting a constant ilow and a high T (such that 10^5/T < 70). As you can see, when I do that convert_u_to_temp fails to converge. So at this point, I decided to give up.

The kludge hack I’ve implemented is in star_maker5.src, if y returns as zero, I set x to zero. This will set pstar as zero, and the condition to make stars, (if a random number is smaller than pstar) will not be met. If y is not zero, I allow the normal progression to happen. I’m not certain if this works or not physically, but it prevents a crash, which is something.

Here is the source of the diagram above, which uses GraphViz:

digraph G {

    subgraph cluster_0 {
        style=filled;
        color=lightgrey;
        sm50 -> sm51 [ label=" because x=NaN"];
        sm51 -> sm52 [ label=" because y=0"];
        sm52 -> sm53 [label=" because coolrate=0"];
        label = "star_maker5.src";
    }

sm53->GSPH1 [label=" because coolingrate=0"];

    subgraph cluster_1 {
        style=filled;
        color=lightgrey;
        GSPH1;
        label = "Grid_StarParticleHandler.C";
    }

GSPH1->GGCR1 [label = " because (Heat-Lambda)=0"];

    subgraph cluster_2 {
        style=filled;
        color=lightgrey;
        GGCR1->GGCR2 [label = " because Lambda=0"];
        GGCR2->GGCR4 [label = " because LambdaFF=0\n and if(logT < Tmax)==TRUE"];
        GGCR2->GGCR3 [label = " because LambdaCmptn=0\n and if(logT < Tmax)\n and if(ComovingCoordinates)==TRUE"];
        label = "Grid_GadgetCoolingRate.C";
    }

"sm50" [
label = "pstar=NaN \n pstar=(bmass/smthresh)*(1.0-exp(-(1.0-beta)*x*dt/tstar)) !Eq(39)"
shape = "record"
];

"sm51" [
label = "x = 1.0 + 1.0/2.0/y - sqrt(1.0/y + 1.0/4.0/y/y)    !Eq(18)"
shape = "record"
];

"sm52" [
label = "y = tstar*coolrate(i,j,k)/d(i,j,k)/dble(d1)/(beta*usn-(1.0-beta)*temp(i,j,k))      !Eq(16)"
shape="record"
];

"sm53" [
label = "subroutine star_maker5(nx, ny, nz, \n d, dm, temp, coolrate, u, v, w, cooltime, \n dt, r, metal, dx, t, z, procnum, \n d1, x1, v1, t1, \n nmax, xstart, ystart, zstart, ibuff, \n imetal, imethod, mintdyn,\n odthresh, masseff, smthresh, level, np,\n xp, yp, zp, up, vp, wp,\n mp, tdp, tcp, metalf)"
shape="record"
];

"GSPH1" [
label = "coolingrate[coolindex] = GadgetCoolingRate\n (log10(temperature[coolindex]), cgsdensity, electronguessptr, zred);"
shape="record"
];

"GGCR1" [
label = "return (Heat - Lambda);"
shape="record"
];

"GGCR2" [
label ="Lambda = LambdaFF + LambdaCmptn;"
shape="record"
];

"GGCR3" [
label = "LambdaCmptn = 5.65e-36 * ne * (T - 2.73 * (1. + redshift)) * POW(1. + redshift, 4.) / nHcgs;"
shape="record"
];

"GGCR4" [
label = "LambdaFF = bff * (nHp + nHep + 4 * nHepp) * ne;"
shape="record"
];

GGCR4->GGfaar1 [label = " because ne=0"];
GGCR3->GGfaar1 [label = " because ne=0"];

    subgraph cluster_3 {
        style=filled;
        color=lightgrey;
        GGfaar1->GGfaar2 [label = " because nHep and nHepp=0" color="red", fontcolor="red"];
        GGfaar1->GGfaar3 [label = " because nHp=0", color="blue", fontcolor="blue"];
        GGfaar3->GGfaar4 [label = " because nH0=1", color="orange", fontcolor="orange"];
        GGfaar4->GGfaar5 [label = " because gJH0ne=0", color="purple", fontcolor="purple"];
        GGfaar4->GGfaar6 [label = " because geH0=0", color="green", fontcolor="green"];
        GGfaar2->GGfaar5 [label = " because gJHe0ne=0", color="turquoise", fontcolor="turquoise"];
        GGfaar2->GGfaar7 [label = " because geHe0=0", color="magenta", fontcolor="magenta"];
        label = "Grid_Gadgetfindabundancesandrates.C";
    }

"GGfaar1" [
label ="ne = nHp + nHep + 2 * nHepp;    /* eqn (38) */"
shape="record"
];

"GGfaar2" [
label = "if((gJHe0ne + geHe0) \<= SMALLNUM) /* no ionization at all */ \n \{\n nHep = 0.0;\n nHepp = 0.0;\n\}"
shape="record"
];

"GGfaar3" [
label = "nHp = 1.0 - nH0;       /* eqn (34) */"
shape="record"
];

"GGfaar4" [
label = "nH0 = aHp / (aHp + geH0 + gJH0ne); /* eqn (33) */"
shape="record"
];

"GGfaar5" [
label = "if(necgs \<= 1.e-25 \|\| J_UV == 0)\n\{gJH0ne = gJHe0ne = gJHepne = 0;\}"
shape="record"
];

"GGfaar6" [
label = "geH0 = flow * GammaeH0[j] + fhi * GammaeH0[j + 1];"
shape="record"
];

"GGfaar7" [
label = "geHe0 = flow * GammaeHe0[j] + fhi * GammaeHe0[j + 1];"
shape="record"
];

GGfaar5->GIPT1 [label=" because J_UV=0"];
    subgraph cluster_4 {
        style=filled;
        color=lightgrey;
        GIPT1;
        GIPT1->GIPT2 [label=" because ilow=nheattab"];
        GIPT2->GIPT3 [label=" because logz \>\>inlog[z]\n(logz = 2 at redshift 99)"];
        label = "GadgetIonizeParamsTable.C";
    }

"GIPT1" [
label ="if(logz \> inlogz[nheattab - 1] \|\| gH0[ilow] == 0 \|\| gH0[ilow + 1] == 0 \|\| nheattab == 0)\n\{\ngJHe0 = gJHep = gJH0 = 0;\nepsHe0 = epsHep = epsH0 = 0;\nJ_UV = 0;\nreturn;\n\}"
shape="record"
];

"GIPT2" [
label="for(i = 0; i \< nheattab; i++)\n\{\nif(inlogz[i] \< logz)\nilow = i;\nelse\nbreak;\}"
shape="record"
];

"GIPT3" [
label ="logz = log10(redshift + 1.0);"
shape="record"
];

GIPT2->TC1 [label = " because 0\<=inlogz\<=0.85\n(first column)"];
GIPT1->TC1 [label = " because gH0[nheattab+1]=0"];

    subgraph cluster_5 {
        style=filled;
        color=lightgrey;
        TC1;
        label = "TREECOOL";
    }

"TC1" [
label = "0.000 3.03516e-14 1.37296e-14 3.04873e-16 1.74434e-25 1.76233e-25 1.00198e-26\n0.005 3.20557e-14 1.47386e-14 3.14717e-16 1.85463e-25 1.87090e-25 1.03701e-26\n0.010 3.37379e-14 1.57232e-14 3.24518e-16 1.96306e-25 1.97710e-25 1.07195e-26\n.\n.\n.\n0.845 5.61822e-16 4.75209e-16 1.16223e-17 4.00000e-27 7.79785e-27 3.51094e-28\n0.850 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00\n .\n.\n.\n(blank, line nheattab)\n (blank, line nheattab+1)"
shape="record"
];

GGfaar6->GMCT1 [label=" because GammaeH0=0\n because 157809.1/T \> 70"];
GGfaar7->GMCT1 [label=" because GammaeHe0=0\n because 285335.4/T \> 70"];

    subgraph cluster_6 {
        style=filled;
        color=lightgrey;
        GMCT1->GMCT2 [label=" because T is too small"];
        GMCT2->GMCT3 [label=" because Tmin is negative\n because MINGASTEMP=0.1"];
        GMCT3->GMCT4 [label=" so GammaeH0 and GammaeHe0 default to zero."];
        label = "GadgetMakeCoolingTable.C";
    }

"GMCT1" [
label ="if(157809.1 / T \< 70)\nGammaeH0[i] = 5.85e-11 * sqrt(T) * exp(-157809.1 / T) * Tfact;\nif(285335.4 / T \< 70)\nGammaeHe0[i] = 2.38e-11 * sqrt(T) * exp(-285335.4 / T) * Tfact;"
shape="record"
];

"GMCT2" [
label ="T = POW(10.0, Tmin + deltaT * i);"
shape="record"
];

"GMCT3" [
label ="  if(MINGASTEMP \> 0.0)\nTmin = log10(0.1 * MINGASTEMP);\nelse\nTmin = 1.0;"
shape="record"
];

"GMCT4" [
label = "AlphaHp[i] = AlphaHep[i] = AlphaHepp[i] = Alphad[i] = GammaeH0[i] = GammaeHe0[i] = GammaeHep[i] = 0;"
shape="record"
];

"TEST1" [
label = "If I force T to be big, and I set a manual ilow, I get \"failed to converge in convert_u_to_temp()\""
shape="record"]

TEST1->GIPT2 [color="green"];
TEST1->GGcutt1 [color="blue"];
TEST1->GMCT2 [color="red"];

    subgraph cluster_7 {
        style=filled;
        color=lightgrey;
        GGcutt1;
        label = "Grid_Gadgetconvertutotemp.C";
    }

"GGcutt1" [
label ="if(iter \>= MAXITER)   /* if it breaks! */\n\{\nprintf(\"failed to converge in convert_u_to_temp()\n\");\nprintf(\"u_input= %e\\nrho_input=%e\\n ne_input=%e\\n\", u_input, rho_input, ne_input);\nprintf(\"DoCool_u_old_input=%e\\nDoCool_rho_input= %e\\nDoCool_dt_input= %e\\nDoCool_ne_guess_input= %e\\n\",\nDoCool_u_old_input, DoCool_rho_input, DoCool_dt_input, DoCool_ne_guess_input);\n\}  /* need to think about how to do error-checking here! */"
shape="record"
];

TEST1->GIVEUP;

"GIVEUP" [
label ="So I give up!"
shape="tripleoctagon"
];

}

I have been getting NaNs at this line in star_maker5.src:

              if(dummy1 .gt. pstar) goto 10

where pstar = NaN, which comes from this line:

c	       Calculate the fraction of mass in cold clouds
               x = 1.0 + 1.0/2.0/y - sqrt(1.0/y + 1.0/4.0/y/y)	!Eq(18)

because y is zero. y is zero because in this, coolrate is zero:

               y = tstar*coolrate(i,j,k)/d(i,j,k)/dble(d1)/
     &              (beta*usn-(1.0-beta)*temp(i,j,k))		!Eq(16)

coolrate is zero because in Grid_StarParticleHandler.C, coolingrate is zero. coolingrate is zero because in Grid_GadgetCoolingRate.C ne is zero.

ne is zero because ne is zero in Gadgetfindabundancesandrates.C. ne is zero because nHP, nHep and nHepp are zero.

nHep and nHepp are zero because gJHe0ne + geHe0 <= SMALLNUM.

nHp is zero because nH0=1. nH0=1 because geH0 and gJHne = 0.

gJHne = 0 because J_UV = 0.

geH0 is zero because GammaeH0[j]=0.

GammaeH0[j] is zero because T is too low GadgetMakeCoolingTable.C. T is small because Tmin is negative. Tmin is negative because MINGASTEMP is 0.1, and Tmin = log10(0.1 * MINGASTEMP), so we get Tmin = -2.

J_UV = 0 because in GadgetIonizeParamsTable.C, this loop for high redshifts always goes to nheattab

  logz = log10(redshift + 1.0);
  ilow = 0;
  for(i = 0; i < nheattab; i++)
    {
      if(inlogz[i] < logz)
	ilow = i;
      else
	break;
    }

but in the same file we have this:

  if(logz > inlogz[nheattab - 1] || gH0[ilow] == 0 || gH0[ilow + 1] == 0 || nheattab == 0)
    {
      gJHe0 = gJHep = gJH0 = 0;
      epsHe0 = epsHep = epsH0 = 0;
      J_UV = 0;
      return;
    }
  else
    J_UV = 1.e-21;		/* irrelevant as long as it's not 0 */

where if ilow = nheattab, gH0[ilow+1] will always be zero because that line in TREECOOL is zero (or equivalently non-existent).

Yesterday I added my galaxy merging code into yt.

Right now I have it more or less working, but here are the things I can think of that need to be fixed:

- I have the last singular halo hardwired to be the zero (largest) halo. I should make it a choice of which final halo to start with, or a set of final haloes to consider (which could be all the haloes).
-I need to see what’s wrong with the topmost (first) graphviz boxes. I think it has to do with the colons and curly braces.

Using a list of galaxy-cluster sized haloes generated using the halo finder VOBOZ (Neyrinck, Gnedin and Hamilton (2005) (MNRAS 356, 1222), I will create mock catalogs based on the method of Yan, Martin and White (2003) (ApJ, 598:848). This method uses a Halo Occupation Distribution (HOD) and Conditional Luminosity Function (CLF) model to ‘place’ galaxies inside the large haloes. By carefully adjusting the HOD and CLF fitting parameters, the mock galaxy statistics can be made to match observed galaxy mass and luminosity functions. A preliminary educated guess is that this process will create several hundred thousand mock galaxies.

Next, using a density-luminosity binning technique described in Gerke et. al. (2007) (MNRAS 376, 1425), I will semi-randomly assign observed galaxies from Sloan or DEEP2 onto each mock galaxy. At redshifts near zero, I will simply use Sloan galaxies, and near a redshift of one, DEEP2 galaxies. For redshifts between Sloan and DEEP2, and past DEEP2, I will use a simple linear fit using the results of Blanton (2006) (ApJ 648:268) to adjust the colors and luminosities of the mock galaxies for redshift. I will output the final mock catalog as a SQLite database which will contain galaxy positions, luminosities, colors and Sloan/DEEP2 object IDs.

I installed yt on my Mac

I’m next going to try to put my galaxy merging stuff into yt.

I’m also hunting down a Gadget bug in enzo.

I need Rick to give me:

For a given tile the list of haloes and particles.

There. It’s written up.