Commits

Matthew Turk  committed d934dbf

We do need the MPI Barrier. For testing runs, I am also projecting Ones. I've
added a join_proj script which can join the on-disk files, although it is
presently serial-only.

  • Participants
  • Parent commits ec2cac3

Comments (0)

Files changed (2)

File join_proj.py

+from yt.mods import *
+from yt.utilities.amr_utils import QuadTree, merge_quadtrees
+import glob
+import h5py
+
+class PFMock(dict):
+    domain_left_edge = na.zeros(3, dtype='float64')
+    domain_right_edge = na.ones(3, dtype='float64')
+pf = PFMock()
+class ProjMock(dict):
+    pass
+
+NN = 32
+cycle = 515
+fields = ["Density", "HII_Density", "HeI_Density", "Metal_Density", "Ones"]
+
+projs = []
+for ax in range(3):
+    qt_w = QuadTree(na.ones(2, "int64") * NN, 5)
+    print "Axis", ax
+    for fn in sorted(glob.glob("cycle%08i/unweighted_*.h5" % cycle)):
+        print fn
+        f = h5py.File(fn)
+        qt_a = QuadTree(na.ones(2, "int64") * NN, 5)
+        qt_a.frombuffer(f["/axis_%02i/refined" % ax][:],
+                        f["/axis_%02i/values" % ax][:],
+                        f["/axis_%02i/weight" % ax][:], "integrate")
+        merge_quadtrees(qt_w, qt_a)
+    # Now we construct our dictionary of values
+    coord_data, field_data, weight_data, dxs = [], [], [], []
+    level = 0
+    while 1:
+        npos, nvals, nwvals = qt_w.get_all_from_level(level)
+        if len(npos) == 0: break
+        coord_data.append(npos)
+        nvals /= nwvals[:,None]
+        field_data.append(nvals)
+        weight_data.append(nwvals)
+        dxs.append(1.0/(NN*2.0**level) * na.ones(nvals.shape[0], dtype="float64"))
+        level += 1
+    coord_data = na.concatenate(coord_data, axis=0).transpose()
+    field_data = na.concatenate(field_data, axis=0).transpose()
+    weight_data = na.concatenate(weight_data, axis=0).transpose()
+    dxs = na.concatenate(dxs, axis=0).transpose()
+    proj = ProjMock()
+    proj.pf = pf
+    proj['px'] = (coord_data[0,:] + 0.5) * dxs
+    proj['py'] = (coord_data[1,:] + 0.5) * dxs
+    proj['pdx'] = dxs / 2.0
+    proj['pdy'] = dxs / 2.0
+    proj['weight_field'] = weight_data
+    for field, arr in zip(fields, field_data):
+        proj[field] = arr
+    projs.append(proj)
+    proj.axis = ax
+    frb = FixedResolutionBuffer(proj, (0.0, 1.0, 0.0, 1.0), (2048, 2048))
+    write_image(na.log10(frb["Density"]), "hi_%02i.png" % ax)

File user_script.py

     ( ["Density", "Temperature", "HII_Fraction",
        "H2I_Fraction", "Metal_Density"],
       "Density"),
-    ( ["Density", "HII_Density", "HeI_Density", "Metal_Density"],
+    ( ["Density", "HII_Density", "HeI_Density", "Metal_Density", "Ones"],
       None),
       )
 
         pf.h
     if root and not os.path.exists("%s" % pf):
         os.makedirs("%s" % pf)
+    MPI.COMM_WORLD.Barrier()
     with time_region("Projections (weighted)"):
         non_reduced_proj(pf, projs[0][0], projs[0][1], "weighted")
     with time_region("Projections (unweighted)"):
         f.create_dataset("axis_%02i/values" % (ax), data=buf[1])
         f.create_dataset("axis_%02i/weight" % (ax), data=buf[2])
         f.flush()
-    f["/"].attrs["fields"] = fields
-    f["/"].attrs["weights"] = str(weight)
-    f["/"].attrs["pf_name"] = str(pf)
-    f["/"].attrs["redshift"] = pf.current_redshift
-    f["/"].attrs["time_seconds"] = pf.current_time * pf["Time"]
+    f.create_group("/attributes")
+    f["/attributes"].attrs["fields"] = fields
+    f["/attributes"].attrs["weights"] = str(weight)
+    f["/attributes"].attrs["pf_name"] = str(pf)
+    f["/attributes"].attrs["redshift"] = pf.current_redshift
+    f["/attributes"].attrs["time_seconds"] = pf.current_time * pf["Time"]
+    f.flush()
     f.close()
     del f