Add a new test.

3 files changed, 67 insertions, 150 deletions

diff --git a/Makefile b/Makefile
index a2eebed..88032bf 100644
--- a/Makefile
+++ b/Makefile
@@ -33,7 +33,7 @@ clean:
 	-rm -f *.o *.d $(TARGET)
 
 test: $(TARGET)
-	LD_LIBRARY_PATH=. ./test.py
+	LD_LIBRARY_PATH=. PYTHONPATH=. ./tests/ham_constraint.py
 
 -include $(OBJS:.o=.d)
 
diff --git a/test.py b/test.py
deleted file mode 100755
index f6e7b51..0000000
--- a/test.py
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/usr/bin/python
-
-import brill_data
-import numpy as np
-import sys
-
-def eval_residual_coeffs(bd, x, z):
-    X, Z = np.meshgrid(x, z)
-
-    psi     = bd.eval_psi(x, z)
-    d2psi_x = bd.eval_psi(x, z, (2, 0))
-    d2psi_z = bd.eval_psi(x, z, (0, 2))
-    dpsi_x  = bd.eval_psi(x, z, (1, 0))
-
-    d2psi_y = np.where(np.abs(X) < 1e-12, d2psi_x, dpsi_x / X)
-
-    d2q = bd.eval_q(x, z, (2, 0)) + bd.eval_q(x, z, (0, 2))
-
-    d2psi = d2psi_x + d2psi_y + d2psi_z
-
-    return d2psi + 0.25 * psi * d2q
-
-coeffs     = (10, 15, 20, 25, 30, 40, 50, 60, 80, 100)
-amplitudes = (1, 2, 4, 6)
-
-x  = np.linspace(0, 10, 100)
-z    = np.linspace(0, 10, 100)
-
-err_coeffs = [
-    [ 1, 0.01, 0.001, 5e-7, 1e-8, 1e-11, 1e-12, 5e-13, 1e-13, 1e-13 ], # amplitude 1
-    [ 5, 0.1,  5e-4,  5e-7, 5e-8, 5e-11, 5e-12, 5e-11, 1e-13, 1e-13 ], # amplitude 2
-    [ 5, 1,    5e-3,  1e-6, 5e-8, 1e-10, 1e-11, 5e-12, 5e-13, 5e-13 ], # amplitude 4
-    [ 5, 15,   1e-2,  1e-5, 1e-7, 5e-10, 5e-11, 5e-12, 5e-13, 1e-12 ], # amplitude 6
-]
-
-sys.stderr.write('Testing convergence of the residual calculated by spectral derivatives\n')
-for amplitude, err_row in zip(amplitudes, err_coeffs):
-    sys.stderr.write('amplitude: %d\n' % amplitude)
-
-    for coeff, err_ref in zip(coeffs, err_row):
-        sys.stderr.write('coeffs: %dx%d ' % (coeff, coeff))
-
-        BD = brill_data.BrillData(nb_coeffs = [coeff, coeff], amplitude = amplitude)
-
-        error = np.max(np.abs(eval_residual_coeffs(BD, x, z)))
-        sys.stderr.write('calculated error %g, reference %g\n' % (error, err_ref))
-        if error > err_ref:
-            sys.stderr.write('Test FAILED\n')
-            sys.exit(1)
-
-
-def eval_residual_fd(bd, x, y, z):
-    Z, Y, X = np.meshgrid(z, y, x, indexing = 'ij')
-
-    psi = np.empty_like(Z)
-    d2q = np.empty_like(Z)
-    for i, yy in enumerate(y):
-        rho = np.sqrt(x ** 2 + yy ** 2)
-        stuff = bd.eval_psi(rho, z)
-        psi[:, i] = bd.eval_psi(rho, z)
-        d2q[:, i] = bd.eval_q(rho, z, (2, 0)) + bd.eval_q(rho, z, (0, 2))
-
-    dpsi_z, dpsi_y, dpsi_x = np.gradient(psi, z[1] - z[0], y[1] - y[0], x[1] - x[0])
-    d2psi_z, foo, foo = np.gradient(dpsi_z, z[1] - z[0], y[1] - y[0], x[1] - x[0])
-    foo, d2psi_y, foo = np.gradient(dpsi_y, z[1] - z[0], y[1] - y[0], x[1] - x[0])
-    foo, foo, d2psi_x = np.gradient(dpsi_x, z[1] - z[0], y[1] - y[0], x[1] - x[0])
-
-    d2psi = d2psi_x + d2psi_y + d2psi_z
-
-    return (d2psi + 0.25 * psi * d2q)[2:-2, 2:-2, 2:-2]
-    return (d2psi + 0.25 * psi * d2q)[2:-2, 2:-2, 2:-2]
-
-points = [ 100,  500,  1000 ]
-err_fd = [ 1e-2, 5e-4, 1e-4 ]
-sys.stderr.write('Testing convergence of the residual calculated by finite difference derivatives\n')
-BD = brill_data.BrillData(nb_coeffs = [50, 50], amplitude = 5)
-for point, err_ref in zip(points, err_fd):
-    x    = np.linspace(-1, 1, point)
-    y    = np.linspace(-3 * (x[1] - x[0]), 3 * (x[1] - x[0]), 7)
-    z    = np.linspace(-1, 1, point)
-    sys.stderr.write('number of points: %d\n' % point)
-
-    error = np.max(np.abs(eval_residual_fd(BD, x, y, z)))
-    sys.stderr.write('calculated error %g, reference %g\n' % (error, err_ref))
-    if error > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
-
-sys.stderr.write('Testing metric evaluation\n')
-
-rho  = np.linspace(0, 1, 100)
-z    = np.linspace(0, 1, 100)
-
-err_metric = [
-    [ 5e-3, 5e-4, 1e-2, 5e-3, 1e-3 ],    # component 00
-    [ 5e-3, 5e-4, 1e-2, 5e-3, 1e-3 ],    # component 11
-    [ 5e-4, 5e-5, 5e-4, 1e-4, 1e-4 ],    # component 22
-]
-
-BD = brill_data.BrillData(nb_coeffs = [40, 40], amplitude = 3)
-for c, err_row in zip(xrange(3), err_metric):
-    sys.stderr.write('component %d%d\n' % (c, c))
-
-    metric  = BD.eval_metric(rho, z, (c, c))
-
-    dmetric_rho      = BD.eval_metric(rho, z, (c, c), (1, 0))
-    dmetric_z        = BD.eval_metric(rho, z, (c, c), (0, 1))
-    d2metric_rho_z   = BD.eval_metric(rho, z, (c, c), (1, 1))
-    d2metric_rho     = BD.eval_metric(rho, z, (c, c), (2, 0))
-    d2metric_z       = BD.eval_metric(rho, z, (c, c), (0, 2))
-
-    emetric_z, emetric_rho       = np.gradient(metric,      z[1] - z[0], rho[1] - rho[0])
-    e2metric_rho_z, e2metric_rho = np.gradient(dmetric_rho, z[1] - z[0], rho[1] - rho[0])
-    e2metric_z, foo              = np.gradient(dmetric_z,   z[1] - z[0], rho[1] - rho[0])
-
-    err = np.max(np.abs((dmetric_rho - emetric_rho)[3:-3, 3:-3]))
-    err_ref = err_row[0]
-    sys.stderr.write('calculated error for d/drho %g, reference %g\n' % (err, err_ref))
-    if err > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
-
-    err = np.max(np.abs((dmetric_z - emetric_z)[3:-3, 3:-3]))
-    err_ref = err_row[1]
-    sys.stderr.write('calculated error for d/dz %g, reference %g\n' % (err, err_ref))
-    if err > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
-
-    err = np.max(np.abs((d2metric_rho - e2metric_rho)[3:-3, 3:-3]))
-    err_ref = err_row[2]
-    sys.stderr.write('calculated error for d2/drho2 %g, reference %g\n' % (err, err_ref))
-    if err > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
-
-    err = np.max(np.abs((d2metric_z - e2metric_z)[3:-3, 3:-3]))
-    err_ref = err_row[3]
-    sys.stderr.write('calculated error for d2/dz2 %g, reference %g\n' % (err, err_ref))
-    if err > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
-
-    err = np.max(np.abs((d2metric_rho_z - e2metric_rho_z)[3:-3, 3:-3]))
-    err_ref = err_row[4]
-    sys.stderr.write('calculated error for d2/drhodz %g, reference %g\n' % (err, err_ref))
-    if err > err_ref:
-        sys.stderr.write('Test FAILED\n')
-        sys.exit(1)
diff --git a/tests/ham_constraint.py b/tests/ham_constraint.py
new file mode 100755
index 0000000..b3e42d3
--- /dev/null
+++ b/tests/ham_constraint.py
@@ -0,0 +1,66 @@
+#!/usr/bin/python
+
+import brill_data
+import numpy as np
+
+import sys
+
+def _invert_matrix(mat):
+    oldshape = mat.shape
+    newshape = oldshape[:2] + (np.product(mat.shape[2:]),)
+
+    mat = mat.reshape(newshape).transpose((2, 0, 1))
+    inv = np.linalg.inv(mat)
+    return inv.transpose((1, 2, 0)).reshape(oldshape)
+
+def ham_calc(bd, rho, z):
+    g   = np.empty((3, 3) + x.shape + z.shape)
+    dg  = np.empty((3,) + g.shape)
+    d2g = np.zeros((3,) + dg.shape)
+
+    for i in range(3):
+        for j in range(3):
+            g[i, j] = bd.eval_metric(rho, z, i * 3 + j)
+
+            dg[0, i, j] = bd.eval_metric(rho, z, i * 3 + j, [1, 0])
+            dg[1, i, j] = bd.eval_metric(rho, z, i * 3 + j, [0, 1])
+            dg[2, i, j] = 0.0
+
+            d2g[0, 0, i, j] = bd.eval_metric(rho, z, i * 3 + j, [2, 0])
+            d2g[0, 1, i, j] = bd.eval_metric(rho, z, i * 3 + j, [1, 1])
+            d2g[1, 0, i, j] = d2g[1, 0, i, j]
+            d2g[1, 1, i, j] = bd.eval_metric(rho, z, i * 3 + j, [0, 2])
+
+            d2g[2, 2, i, j] = 0.0
+            d2g[0, 2, i, j] = 0.0
+            d2g[2, 0, i, j] = 0.0
+            d2g[1, 2, i, j] = 0.0
+            d2g[2, 1, i, j] = 0.0
+
+    gu = _invert_matrix(g)
+
+    Gl = 0.5 * (np.einsum('cab...->abc...', dg) + np.einsum('bac...->abc...', dg) - dg)
+    G  = np.einsum('ij...,jkl...->ikl...', gu, Gl)
+    R = (0.5 * (np.einsum('adbc...->abcd...', d2g) + np.einsum('bcad...->abcd...', d2g) - np.einsum('bdac...->abcd...', d2g) - np.einsum('acbd...->abcd...', d2g)) +
+         np.einsum('ead...,ebc...->abcd...', Gl, G) - np.einsum('eac...,ebd...->abcd...', Gl, G))
+
+    R = np.einsum('ijkl...,ik...->jl...', R, gu)
+    R = np.einsum('ij...,ij...', R, gu)
+    return R
+
+tol = 5e-13
+
+x = np.linspace(0.1, 256, 50)
+z = np.linspace(0.1, 256, 50)
+
+coeffs = (80, 30)
+amplitudes = np.linspace(1, 12, 11)
+hmax = np.empty_like(amplitudes)
+for i, a in enumerate(amplitudes):
+    bd = brill_data.BrillData(amplitude = a, nb_coeffs = coeffs)
+    hmax[i] = np.max(np.abs(ham_calc(bd, x, z)))
+
+if np.any(hmax > tol):
+    sys.stderr.write('Large constraint violation: amplitudes %s, H %s\n'
+                     % (amplitudes, hmax))
+    sys.exit(1)