Merge pull request #7 from daviddritschel/write-h5py

Convert script to hdf5 file

scripts/dv

@@ -24,9 +24,9 @@ else
 endif
 
 # Get resolution
-set nx = `grep "nx=" src/parameters.f90 | awk -F= '{print $(NF-1)}' | awk -F, '{print $1}'`
+set nx = `grep "nx=" PS_SRC_DIR/parameters.f90 | awk -F= '{print $(NF-1)}' | awk -F, '{print $1}'`
 
-set ny = `grep "ny=" src/parameters.f90 | awk -F= '{print $(NF)}' | awk -F, '{print $1}'`
+set ny = `grep "ny=" PS_SRC_DIR/parameters.f90 | awk -F= '{print $(NF)}' | awk -F, '{print $1}'`
 
 @ nyp1 = ( $ny + 1 )
 

scripts/image.py

@@ -15,6 +15,7 @@
 import matplotlib as mpl
 from matplotlib import rcParams
 from matplotlib import rc
+from ps_config import nx, ny
 rcParams.update({'figure.autolayout': True})
 warnings.simplefilter("ignore",DeprecationWarning)
 
@@ -98,16 +99,7 @@ def contint(fmin,fmax):
 t = float(t_in or 900.0)
 
 #-----------------------------------------------------------------
-# Work out grid resolution (ng) by reading it from parameters.f90:
-in_file=open('src/parameters.f90','r')
-fread=in_file.readlines()
-for line in fread:
-   if ':: nx=' in line:
-      pline=line
-
-line=pline.split("=")[1]
-nx=int(line.split(",")[0])
-ny=int(pline.split("=")[2])+1
+ny = ny + 1
 
 #-----------------------------------------------------------------
 # Open ene.asc file in one directory to get time between frames:

scripts/robert_image.py

@@ -15,6 +15,7 @@
 import matplotlib as mpl
 from matplotlib import rcParams
 from matplotlib import rc
+from ps_config import nx, ny
 rcParams.update({'figure.autolayout': True})
 warnings.simplefilter("ignore",DeprecationWarning)
 
@@ -98,16 +99,7 @@ def contint(fmin,fmax):
 t = float(t_in or 900.0)
 
 #-----------------------------------------------------------------
-# Work out grid resolution (ng) by reading it from parameters.f90:
-in_file=open('src/parameters.f90','r')
-fread=in_file.readlines()
-for line in fread:
-   if ':: nx=' in line:
-      pline=line
-
-line=pline.split("=")[1]
-nx=int(line.split(",")[0])
-ny=int(pline.split("=")[2])+1
+ny = ny + 1
 
 #-----------------------------------------------------------------
 # Open ene.asc file in one directory to get time between frames:

scripts/straka_image.py

@@ -15,6 +15,7 @@
 import matplotlib as mpl
 from matplotlib import rcParams
 from matplotlib import rc
+from ps_config import nx, ny
 rcParams.update({'figure.autolayout': True})
 warnings.simplefilter("ignore",DeprecationWarning)
 
@@ -97,17 +98,7 @@ def contint(fmin,fmax):
 t_in = input(' Time to show (default 900)? ')
 t = float(t_in or 900.0)
 
-#-----------------------------------------------------------------
-# Work out grid resolution (ng) by reading it from parameters.f90:
-in_file=open('src/parameters.f90','r')
-fread=in_file.readlines()
-for line in fread:
-   if ':: nx=' in line:
-      pline=line
-
-line=pline.split("=")[1]
-nx=int(line.split(",")[0])
-ny=int(pline.split("=")[2])+1
+ny = ny +1
 
 #-----------------------------------------------------------------
 # Open ene.asc file in one directory to get time between frames:

scripts/write_h5.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+import h5py
+import numpy as np
+import os
+from ps_config import nx, ny, ellx, ymin, ymax
+
+print("Parameters from parameters.f90:")
+print("nx   = ", nx)
+print("ny   = ", ny)
+print("ellx = ", ellx)
+print("ymin = ", ymin)
+print("ymax = ", ymax)
+
+#-----------------------------------------------------------------
+# Open ene.asc file in one directory to get time between frames:
+in_file=open('evolution/ecomp.asc','r')
+time, ekin, epot, etot = np.loadtxt(in_file,dtype=float,unpack=True)
+in_file.close()
+nt = len(time)
+
+dset = {
+    'bb': 'buoyancy',
+    'zz': 'vorticity'
+}
+
+N = nx * (ny + 1)
+
+h5file = h5py.File('ps_fields.hdf5', 'w')
+
+h5file.attrs['nsteps'] = nt
+
+box = h5file.create_group('box')
+box.attrs['extent'] = (ellx, ymax - ymin)
+box.attrs['origin'] = (-ellx * 0.5, ymin)
+box.attrs['ncells'] = (np.int32(nx), np.int32(ny))
+
+for frame in range(nt):
+    group = h5file.create_group('step#' + str(frame).zfill(10))
+
+    for field in ['bb', 'zz']:
+        fname = os.path.join('evolution', field + '.r4')
+
+        in_file = open(fname,'r')
+        raw_array = np.fromfile(in_file,dtype=np.float32)
+        in_file.close()
+
+        Z = np.empty([nx, ny+1])
+        Z = raw_array[frame*(N+1)+1:(frame+1)*(N+1)].reshape(nx, ny+1)
+
+        group[dset[field]] = np.float64(Z.copy())
+
+h5file.close()

src/makefile

@@ -49,6 +49,16 @@ setup:
 install:
 	mkdir -p $(installdir)
 	find -type f -perm /a+x -exec /bin/mv {} $(installdir)/ \;
+	# extract parameters.f90 values
+	cat parameters.f90 | grep -oE nx=[[:digit:]]* > ps_config.py
+	cat parameters.f90 | grep -oE ny=[[:digit:]]* >> ps_config.py
+	cat parameters.f90 | grep -oE ellx=[[:digit:]]*[.]*[[:digit:]]* >> ps_config.py
+	cat parameters.f90 | grep -oE ymin=[[:digit:]]*[.]*[[:digit:]]* >> ps_config.py
+	cat parameters.f90 | grep -oE ymax=[[:digit:]]*[.]*[[:digit:]]* >> ps_config.py
+	mv ps_config.py $(installdir)
+	cp -r ../scripts/* $(installdir)
+	sed -i "s:PS_SRC_DIR:$(shell pwd):g" $(installdir)/dv
+
 
  # Include intialisation and post-processing rules
  # if these exist: