1. Matthew Turk
  2. axisem

Commits

Matthew Turk  committed 2811502

Refactoring a bit

  • Participants
  • Parent commits 833137d
  • Branches default

Comments (0)

Files changed (3)

File axisem_read.py

View file
  • Ignore whitespace
+from yt.mods import *
+import os
+import numpy as np
+import pyximport; pyximport.install()
+import grid_interp
+
+def read_coords(nproc = 4):
+    coords = np.concatenate([np.array([[float(i) for i in line.split()]
+                            for line in open("glob_grid_%04i.dat" % g)])
+                            for g in range(nproc)])
+    print "Creating QT", coords.min(axis=0), coords.max(axis=0)
+    QT = grid_interp.QTree(coords.min(axis=0), coords.max(axis=0))
+    QT.add(coords)
+    return coords, QT
+
+def read_snapshot(coords, QT, sn, nproc=4, dims = (256, 256, 256), ncol = 5):
+    # Now we have our x, y coords in the plane of the simulation.  We need to
+    # interpolate into 3D.
+
+    # It's supposed to be a sphere, so these will define our left edge and right
+    # edge for our grids.
+    mi = coords.min()
+    ma = coords.max()
+
+    # Now we set up our grids.  But only if it's not already on-disk.
+
+    d = np.zeros(dims)
+    arr = np.array(dims, dtype="int64")
+    data = {}
+    print "Allocating memory"
+    for col in range(ncol):
+        data["Column%02i" % col] = np.zeros(dims, dtype="float64")
+    pf = load_uniform_grid(data, arr, 100.0,
+                           bbox = np.array( [[mi, ma], [mi, ma], [mi, ma]] ),
+                           nprocs = 32)
+    del data
+    #pf.h.print_stats()
+    fvals = []
+    print "Reading from disk"
+    for g in range(nproc):
+        arr = np.fromfile("snap_%04i_%04i.dat" % (g, sn), dtype="f")
+        fvals.append(arr.reshape(arr.size / ncol, ncol, order="C"))
+    fvals = np.concatenate(fvals).astype("float64")
+    np.abs(fvals, fvals)
+
+    pb = get_pbar("Interpolating ", len(pf.h.grids))
+    for i, g in enumerate(pf.h.grids):
+        pb.update(i+1)
+        for col in range(ncol):
+            d = pf.h.io.fields[g.id]["Column%02i" % col]
+            QT.interp(g.LeftEdge, g.RightEdge, d, fvals[:,col], ma)
+            #np.clip(g["Column%02i" % col], 1e-20, 1e20,
+            #        g["Column%02i" % col])
+    pb.finish()
+
+    return pf, fvals

File convert.py

View file
  • Ignore whitespace
 import numpy as np
 import pyximport; pyximport.install()
 import grid_interp
+from axisem_read import read_coords, read_snapshot
 
-coords = np.concatenate([np.array([[float(i) for i in line.split()]
-                        for line in open("glob_grid_%04i.dat" % g)]) for g in range(4)])
+for snap in parallel_objects(range(1, 360)):
+    coords, QT = read_coords()
+    pf, fvals = read_snapshot(coords, QT, 200, dims=(512,512,512))
 
-QT = grid_interp.QTree(coords.min(axis=0), coords.max(axis=0))
-QT.add(coords)
-level_counts = QT.stats()
-for i, c in enumerate(level_counts):
-    if c == 0: break
-    print "% 3i % 8i" % (i, c)
+    for col in range(5):
+        proj = pf.h.proj(0, "Column%02i" % col, style="mip")
+        pw = proj.to_pw()
+        pw.save("SNAP_%03i" % snap)
+        pw.set_log("Column%02i" % col, True)
+        pw.save("SNAP_%03i_log" % snap)
 
-# Now we have our x, y coords in the plane of the simulation.  We need to
-# interpolate into 3D.
+#tf = ColorTransferFunction((-8, -3))
+#tf.add_layers(8)
 
-# It's supposed to be a sphere, so these will define our left edge and right
-# edge for our grids.
-mi = coords.min()
-ma = coords.max()
+#c = np.zeros(3, "float64")
+#L = -np.ones(3, "float64")
+#w = pf.domain_width.copy()
 
-# Now we set up our grids.  But only if it's not already on-disk.
-
-d = np.zeros((512,512,512))
-arr = np.array((512,512,512), dtype="int64")
-pf = load_uniform_grid({"Perturbation": d}, arr, 100.0,
-                       bbox = np.array( [[mi, ma], [mi, ma], [mi, ma]] ),
-                       nprocs = 32)
-#pf.h.print_stats()
-values = (coords**2.0).sum(axis=1)**0.5
-
-#pb = get_pbar("Interpolating ", len(pf.h.grids))
-for i, g in enumerate(pf.h.grids):
-    print i
-    #pb.update(i+1)
-    #print g.LeftEdge, g.RightEdge
-    QT.interp(g.LeftEdge, g.RightEdge, g["Perturbation"], values,
-              ma)
-#pb.finish()
-p = SlicePlot(pf, "x", "Perturbation")
-p.save()
+#cam = pf.h.camera(c, L, w, 1024, tf,
+#        fields = ["Column02"], log_fields=[True])
+#s = cam.snapshot("hi.png", 4.0)
+#print s.min(), s.max()
+#nn = np.isnan(s[:,:,0])
+#s[nn,:] = 0.0
+#write_bitmap(s, "hi.png")

File grid_interp.pyx

View file
  • Ignore whitespace
                         if r > max_r: continue
                     if node.n == 0:
                         # We have nothing to interpolate
+                        #print "SKIPPING NODE", r
                         continue
                     #print pos[0], pos[1], pos[2], cpos[0], cpos[1],
-                    #print node.n
+                    #print node.n,
                     #continue
                     ind = self.find_nearest(cpos)
+                    #print ind, values[ind], i, j, k
                     dest[i,j,k] = values[ind]
         #raise RuntimeError