Switch to nr_analysis_axi interpolation.HEADmaster

It is more accurate.

Also, update for python3.

1 files changed, 33 insertions, 17 deletions

diff --git a/revive.py b/revive.py
index dbc3658..67010bf 100755
--- a/revive.py
+++ b/revive.py
@@ -1,4 +1,7 @@
-#!/usr/bin/python
+#!/usr/bin/python3
+
+# This script is intended to be invoked by the ReviveData cactus thorn [1].
+# [1] https://git.khirnov.net/thorns/ReviveData.git/
 
 import argparse
 import sys
@@ -7,6 +10,7 @@ import numpy as np
 from scipy import interpolate
 
 from cactus_utils import simdata
+from nr_analysis_axi import interp
 
 metric_str = {
     'gxx' : '11', 'gyy' : '22', 'gzz' : '33',
@@ -20,6 +24,8 @@ shift_str = {
     'betax' : 'beta1', 'betay' : 'beta2', 'betaz' : 'beta3',
 }
 
+interp_order = 8
+
 parser = argparse.ArgumentParser(description = 'Revive cactus simulation data')
 
 parser.add_argument('simdata_path')
@@ -35,6 +41,11 @@ parser.add_argument('gridsize_z', type = int)
 
 args = parser.parse_args(sys.argv[1:])
 
+x = np.linspace(args.x0, args.x0 + args.dx * (args.gridsize_x - 1), args.gridsize_x)
+z = np.linspace(args.z0, args.z0 + args.dz * (args.gridsize_z - 1), args.gridsize_z)
+
+X, Z = np.meshgrid(x, z)
+
 sd = simdata.SimulationData(args.simdata_path)
 
 if args.load_iteration < 0:
@@ -51,8 +62,7 @@ load_rl = None
 for rl in sd.rl[::-1]:
     s = rl.slice(load_iteration)
     if (s.x[0] <= args.x0 and s.z[0] <= args.z0 and
-        s.x[-1] >= (args.x0 + args.dx * (args.gridsize_x - 1)) and
-        s.z[-1] >= (args.z0 + args.dz * (args.gridsize_z - 1))):
+        s.x[-1] >= x[-1] and s.z[-1] >= z[-1]):
         load_rl = rl
         break
 
@@ -61,21 +71,27 @@ if load_rl is None:
     sys.exit(1)
 
 s = load_rl.slice(load_iteration)
-if varname in metric_str:
-    src = s['gt' + metric_str[varname]][:] / (s['phi'][:] ** 2)
-elif varname in curv_str:
-    src = (s['At' + curv_str[varname]][:] + (1.0 / 3.0) * s['gt' + curv_str[varname]][:] * s['trK'][:]) / (s['phi'][:] ** 2)
-elif varname in shift_str:
-    src = s[shift_str[varname]][:]
-elif varname == 'alp':
-    src = s['alpha'][:]
+src = None
+try:
+    if varname in metric_str:
+        src = s['gt' + metric_str[varname]][:] / (s['phi'][:] ** 2)
+    elif varname in curv_str:
+        src = (s['At' + curv_str[varname]][:] + (1.0 / 3.0) * s['gt' + curv_str[varname]][:] * s['trK'][:]) / (s['phi'][:] ** 2)
+    elif varname in shift_str:
+        src = s[shift_str[varname]][:]
+    elif varname == 'alp' or varname == 'alpha':
+        src = s['alpha'][:]
+    else:
+        sys.stderr.write('Unknown gridfunction: %s\n' % varname)
+        sys.exit(1)
+except IndexError:
+    sys.stderr.write('Gridfunction "%s" not in data, assuming zero\n' % varname)
+    np.zeros(z.shape + x.shape).tofile(sys.stdout.buffer.raw)
+    sys.exit(0)
 
 src = src.squeeze()
 
-src_interp = interpolate.RectBivariateSpline(s.z, s.x, src)
-
-x = np.linspace(args.x0, args.x0 + args.dx * (args.gridsize_x - 1), args.gridsize_x)
-z = np.linspace(args.z0, args.z0 + args.dz * (args.gridsize_z - 1), args.gridsize_z)
+src_interp = interp.Interp2D_C([s.z[0], s.x[0]], [s.z[1] - s.z[0], s.x[1] - s.x[0]], src, interp_order)
 
-res = src_interp(z, x)
-res.tofile(sys.stdout)
+res = src_interp(Z.flatten(), X.flatten())
+res.tofile(sys.stdout.buffer.raw)