From bd37d6d1c90dbfb5ae882c8eb9bb779bd00962c3 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 25 Jul 2017 15:01:29 +0200 Subject: Add a new test. --- Makefile | 2 +- test.py | 149 ------------------------------------------------ tests/ham_constraint.py | 66 +++++++++++++++++++++ 3 files changed, 67 insertions(+), 150 deletions(-) delete mode 100755 test.py create mode 100755 tests/ham_constraint.py 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) -- cgit v1.2.3