From ae210faae56af10e5897fa413e731921d718b8c7 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 2 Dec 2014 19:50:11 +0100 Subject: Add tests. --- Makefile | 5 ++- test.py | 149 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+), 1 deletion(-) create mode 100755 test.py diff --git a/Makefile b/Makefile index 75de828..aec6c0b 100644 --- a/Makefile +++ b/Makefile @@ -13,4 +13,7 @@ $(TARGET): $(OBJECTS) clean: -rm $(OBJECTS) $(TARGET) -.PHONY: clean +test: $(TARGET) + LD_LIBRARY_PATH=. ./test.py + +.PHONY: clean test diff --git a/test.py b/test.py new file mode 100755 index 0000000..5af80bc --- /dev/null +++ b/test.py @@ -0,0 +1,149 @@ +#!/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_rho = 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_rho = 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_rho = 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) -- cgit v1.2.3