add function to compute many point centered differences

bfps/tools.py

@@ -69,4 +69,38 @@ def padd_with_zeros(
     b[n0-m0/2:    , n1-m1/2:    , :m2/2+1] = a[m0-m0/2:    , m1-m1/2:    , :m2/2+1]
     return b
 
+try:
+    import sympy as sp
+
+    def get_fornberg_coeffs(
+            x = None,
+            a = None):
+        N = len(a) - 1
+        d = []
+        for m in range(N+1):
+            d.append([])
+            for n in range(N+1):
+                d[m].append([])
+                for j in range(N+1):
+                    d[m][n].append(sp.Rational(0))
+        d[0][0][0] = sp.Rational(1)
+        c1 = sp.Rational(1)
+        for n in range(1, N+1):
+            c2 = sp.Rational(1)
+            for j in range(n):
+                c3 = a[n] - a[j]
+                c2 = c2*c3
+                for m in range(n+1):
+                    d[m][n][j] = ((a[n] - x)*d[m][n-1][j] - m*d[m-1][n-1][j]) / c3
+            for m in range(n+1):
+                d[m][n][n] = (c1 / c2)*(m*d[m-1][n-1][n-1] - (a[n-1] - x)*d[m][n-1][n-1])
+            c1 = c2
+        coeffs = []
+        for m in range(len(d)):
+            coeffs.append([])
+            for j in range(len(d)):
+                coeffs[-1].append(d[m][N][j])
+        return np.array(coeffs).astype(np.float)
+except ImportError:
+    print('didn\'t find sympy.')
 

tests/test_base.py

@@ -128,30 +128,32 @@ def launch(
               njobs = opt.njobs)
     return c
 
-def acceleration_test(c):
+
+def acceleration_test(c, n = 3):
     import numpy as np
     import matplotlib.pyplot as plt
+    from bfps.tools import get_fornberg_coeffs
     d = c.get_data_file()
     pos = d['particles/tracers4/state'].value
     vel = d['particles/tracers4/velocity'].value
     acc = d['particles/tracers4/acceleration'].value
 
-    num_acc1 = (- vel[ :-2] + vel[2:])/(2*d['parameters/dt'].value*d['parameters/niter_part'].value)
-    num_acc2 = (pos[ :-2] - 2*pos[1:-1] + pos[2:])/((d['parameters/dt'].value*d['parameters/niter_part'].value)**2)
-    num_acc3 = (-vel[4:] + 8*vel[3:-1] - 8*vel[1:-3] + vel[:-4])/(12*d['parameters/dt'].value*d['parameters/niter_part'].value)
-    num_acc4 = (-pos[4:] + 16*pos[3:-1] - 30*pos[2:-2] + 16*pos[1:-3] - pos[:-4])/(12*(d['parameters/dt'].value*d['parameters/niter_part'].value)**2)
-    num_vel = (- pos[ :-2] + pos[2:])/(2*d['parameters/dt'].value*d['parameters/niter_part'].value)
+    fc = get_fornberg_coeffs(0, range(-n, n+1))
+    dt = d['parameters/dt'].value*d['parameters/niter_part'].value
+
+    print [len(range(n-i, vel.shape[0]-i-n)) for i in range(-n, n+1)]
+    num_acc1 = sum(fc[1, n-i]*vel[n-i:vel.shape[0]-i-n] for i in range(-n, n+1)) / dt
+    num_acc2 = sum(fc[2, n-i]*pos[n-i:pos.shape[0]-i-n] for i in range(-n, n+1)) / dt**2
 
-    pid = np.unravel_index(np.argmax(np.abs(num_acc3 - acc[2:-2])), dims = num_acc3.shape)[1]
+    pid = np.unravel_index(np.argmax(np.abs(num_acc1 - acc[n:-n])), dims = num_acc1.shape)[1]
     fig = plt.figure()
     a = fig.add_subplot(111)
     col = ['red', 'green', 'blue']
     for cc in range(3):
-        a.plot(num_acc1[2:, pid, cc], color = col[cc], dashes = (3, 4))
-        a.plot(num_acc2[2:, pid, cc], color = col[cc], dashes = (2, 2))
-        a.plot(num_acc3[1:, pid, cc], color = col[cc])
-        a.plot(num_acc4[1:, pid, cc], color = col[cc], dashes = (3, 5))
-        a.plot(acc[3:, pid, cc], color = col[cc], dashes = (1, 1))
+        a.plot(num_acc1[1:, pid, cc], color = col[cc])
+        a.plot(num_acc2[1:, pid, cc], color = col[cc], dashes = (2, 2))
+        a.plot(acc[n+1:, pid, cc], color = col[cc], dashes = (1, 1))
+    fig.tight_layout()
     fig.savefig(os.path.join(c.work_dir, 'acc_test_{0}.pdf'.format(c.simname)))
     return None