Overhaul of tests sub-package.

tests/analytical/hertzian_dipole_fs/hertzian_dipole_fs.in

@@ -1,8 +0,0 @@
-#title: Hertzian dipole in free-space
-#domain: 0.100 0.100 0.100
-#dx_dy_dz: 0.001 0.001 0.001
-#time_window: 3e-9
-
-#waveform: gaussiandot 1 1e9 mypulse
-#hertzian_dipole: z 0.050 0.050 0.050 mypulse
-#rx: 0.070 0.070 0.070 output Ex Ey Ez Hx Hy Hz

tests/models_basic/2D_ExHyHz/2D_ExHyHz.in

@@ -0,0 +1,8 @@
+#title: 2D test Ex, Hy, Hz components
+#domain: 0.001 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: x 0 0.050 0.050 myWave
+#rx: 0 0.070 0.070

tests/models_basic/2D_EyHxHz/2D_EyHxHz.in

@@ -0,0 +1,8 @@
+#title: 2D test Ey, Hx, Hz components
+#domain: 0.100 0.001 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: y 0.050 0 0.050 myWave
+#rx: 0.070 0 0.070

tests/models_basic/2D_EzHxHy/2D_EzHxHy.in

@@ -0,0 +1,8 @@
+#title: 2D test Ez, Hx, Hy components
+#domain: 0.100 0.100 0.001
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0 myWave
+#rx: 0.070 0.070 0

tests/models_basic/cylinder_Ascan_2D/cylinder_Ascan_2D.in

@@ -0,0 +1,13 @@
+#title: A-scan from a metal cylinder buried in a dielectric half-space
+#domain: 0.240 0.210 0.002
+#dx_dy_dz: 0.002 0.002 0.002
+#time_window: 3e-9
+
+#material: 6 0 1 0 half_space
+
+#waveform: ricker 1 1.5e9 my_ricker
+#hertzian_dipole: z 0.100 0.170 0 my_ricker
+#rx: 0.140 0.170 0
+
+#box: 0 0 0 0.240 0.170 0.002 half_space
+#cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec

tests/models_basic/hertzian_dipole_dispersive/hertzian_dipole_dispersive.in

@@ -0,0 +1,12 @@
+#title: Hertzian dipole in water
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+#rx: 0.070 0.070 0.070
+
+#material: 4.9 0 1 0 myWater
+#add_dispersion_debye: 1 75.2 9.231e-12 myWater
+#box: 0 0 0 0.100 0.100 0.100 myWater

tests/models_basic/hertzian_dipole_fs/hertzian_dipole_fs.in

@@ -0,0 +1,8 @@
+#title: Hertzian dipole in free-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+#rx: 0.070 0.070 0.070

tests/models_basic/hertzian_dipole_fs_analytical/hertzian_dipole_fs_analytical.in

@@ -0,0 +1,8 @@
+#title: Hertzian dipole in free-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+#rx: 0.070 0.070 0.070

tests/models_basic/hertzian_dipole_hs/hertzian_dipole_hs.in

@@ -0,0 +1,11 @@
+#title: Hertzian dipole over a half-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+#rx: 0.070 0.070 0.070
+
+#material: 8 0 1 0 half_space
+#box: 0 0 0 0.100 0.100 0.050 half_space

tests/test_compare_analytical.py

@@ -1,129 +0,0 @@
-# Copyright (C) 2015-2017: The University of Edinburgh
-#                 Authors: Craig Warren and Antonis Giannopoulos
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
-import argparse
-import os
-import sys
-
-import h5py
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-
-from tests.analytical_solutions import hertzian_dipole_fs
-
-"""Plots a comparison of analytical solutions and given simulated output."""
-
-# Parse command line arguments
-parser = argparse.ArgumentParser(description='Plots a comparison of analytical solutions and given simulated output.', usage='cd gprMax; python -m tests.test_compare_analytical modelfile')
-parser.add_argument('modelfile', help='name of model output file including path')
-args = parser.parse_args()
-
-# Model results
-f = h5py.File(args.modelfile, 'r')
-path = '/rxs/rx1/'
-availablecomponents = list(f[path].keys())
-
-floattype = f[path + availablecomponents[0]].dtype
-iterations = f.attrs['Iterations']
-dt = f.attrs['dt']
-dxdydz = f.attrs['dx, dy, dz']
-time = np.linspace(0, 1, iterations)
-time *= (iterations * dt)
-rxpos = f[path].attrs['Position']
-txpos = f['/srcs/src1/'].attrs['Position']
-rxposrelative = ((rxpos[0] - txpos[0]), (rxpos[1] - txpos[1]), (rxpos[2] - txpos[2]))
-model = np.zeros((iterations, len(availablecomponents)), dtype=floattype)
-
-# Analytical solution of a dipole in free space
-analytical = hertzian_dipole_fs(iterations, dt, dxdydz, rxposrelative)
-
-# Read modelled fields and calculate differences
-threshold = 1e-4 # Threshold, below which ignore differences
-diffs = np.zeros((iterations, len(availablecomponents)), dtype=floattype)
-for index in range(len(availablecomponents)):
-    model[:,index] = f[path + availablecomponents[index]][:]
-    max = np.amax(np.abs(analytical[:,index]))
-    if max < threshold:
-        diffs[:,index] = 0
-        diffsum = 0
-        print('Detected differences of less than threshold {}, when comparing {} field component, therefore set as zero.'.format(threshold, availablecomponents[index]))
-    else:
-        diffs[:,index] = (np.abs(analytical[:,index] - model[:,index]) / max) * 100
-        diffsum = (np.sum(np.abs(analytical[:,index] - model[:,index])) / np.sum(np.abs(analytical[:,index]))) * 100
-    print('Total differences in field component {}: {:.1f}%'.format(availablecomponents[index], diffsum))
-
-f.close()
-
-# Plot modelled and analytical solutions
-fig1, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num=args.modelfile + ' versus analytical solution', figsize=(20, 10), facecolor='w', edgecolor='w')
-gs1 = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3)
-
-for index in range(len(availablecomponents)):
-    i = int(index % 3)
-    j = int((index - i) / 3 % 2)
-    ax = plt.subplot(gs1[i, j])
-    line1, = ax.plot(time, model[:,index],'r', lw=2, label='Model')
-    line2, = ax.plot(time, analytical[:,index],'r', lw=2, ls='--', label='Analytical')
-    ax.set_ylim(1.1 * np.amin(np.amin(model[:, 0:3], axis=1)), 1.1 * np.amax(np.amax(model[:, 0:3], axis=1)))
-    
-    if index > 2:
-        plt.setp(line1, color='g')
-        plt.setp(line2, color='g')
-        ax.set_ylim(1.1 * np.amin(np.amin(model[:, 3:6], axis=1)), 1.1 * np.amax(np.amax(model[:, 3:6], axis=1)))
-
-    ax.set_xlim(0, time[-1])
-    ax.grid()
-    ax.legend()
-
-# Set axes labels, limits and turn on grid
-ylabels = ['Ex, field strength [V/m]', 'Ey, field strength [V/m]', 'Ez, field strength [V/m]', 'Hx, field strength [A/m]', 'Hy, field strength [A/m]', 'Hz, field strength [A/m]']
-[ax.set_ylabel(ylabels[index]) for index, ax in enumerate(fig1.axes)]
-
-# Plot differences of modelled and analytical solutions
-fig2, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), num=args.modelfile + ' versus analytical solution differences', figsize=(20, 10), facecolor='w', edgecolor='w')
-gs2 = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3)
-
-for index in range(len(availablecomponents)):
-    i = int(index % 3)
-    j = int((index - i) / 3 % 2)
-    ax = plt.subplot(gs2[i, j])
-    line1, = ax.plot(time, diffs[:, index],'r', lw=2)
-    ax.set_ylim(0, 1.1 * np.amax(np.amax(diffs[:, 0:3], axis=1)))
-    
-    if index > 2:
-        plt.setp(line1, color='g')
-        ax.set_ylim(0, 1.1 * np.amax(np.amax(diffs[:, 3:6], axis=1)))
-
-    ax.set_ylim(0, 2)
-    ax.set_xlim(0, time[-1])
-    ax.grid()
-
-# Set axes labels, limits and turn on grid
-ylabels = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz']
-ylabels = [ylabel + ', percentage difference [%]' for ylabel in ylabels]
-[ax.set_ylabel(ylabels[index]) for index, ax in enumerate(fig2.axes)]
-
-# Save a PDF/PNG of the figure
-savename = os.path.abspath(os.path.dirname(args.modelfile)) + os.sep + os.path.splitext(os.path.split(args.modelfile)[1])[0] + '_vs_analytical'
-#fig1.savefig(savename + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-#fig2.savefig(savename + '_diffs.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-fig1.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
-fig2.savefig(savename + '_diffs.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
-
-plt.show()

tests/test_compare_numerical.py

@@ -1,128 +0,0 @@
-# Copyright (C) 2015-2017: The University of Edinburgh
-#                 Authors: Craig Warren and Antonis Giannopoulos
-#
-# This file is part of gprMax.
-#
-# gprMax is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# gprMax is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
-
-import os
-import sys
-
-import h5py
-import numpy as np
-import matplotlib.pyplot as plt
-
-"""Compare field outputs
-    
-Usage:
-    cd gprMax
-    python -m tests.test_compare_numerical path_to_file1 path_to_file2
-    
-"""
-
-filename1 = sys.argv[1]
-filename2 = sys.argv[2]
-path = '/rxs/rx1/'
-# Key refers to subplot location
-fields = {0: 'Ex', 2: 'Ey', 4: 'Ez', 1: 'Hx', 3: 'Hy', 5: 'Hz'}
-plotorder = list(fields.keys())
-
-# File 1 results
-f = h5py.File(filename1, 'r')
-floattype = f[path + 'Ex'].dtype
-data1 = np.zeros((f.attrs['Iterations'], 6), dtype=floattype)
-time1 = np.zeros((f.attrs['Iterations']), dtype=floattype)
-time1 = np.arange(0, f.attrs['dt'] * f.attrs['Iterations'], f.attrs['dt']) / 1e-9
-for ID, name in fields.items():
-    data1[:,ID] = f[path + str(name)][:]
-f.close()
-
-# File 2 results
-f = h5py.File(filename2, 'r')
-data2 = np.zeros((f.attrs['Iterations'], 6), dtype=floattype)
-time2 = np.zeros((f.attrs['Iterations']), dtype=floattype)
-time2 = np.arange(0, f.attrs['dt'] * f.attrs['Iterations'], f.attrs['dt']) / 1e-9
-for ID, name in fields.items():
-    data2[:,ID] = f[path + str(name)][:]
-f.close()
-
-# Differences
-# In case there is any difference in the number of iterations, take the smaller
-timesmallest = np.amin((time2.shape, time1.shape))
-fieldssmallest = np.amin((data2.shape[0], data1.shape[0]))
-
-threshold = 1e-4 # Threshold, below which ignore differences
-diffs = np.zeros((fieldssmallest, 6), dtype=floattype)
-for ID, name in fields.items():
-    max = np.amax(np.abs(data1[:fieldssmallest,ID]))
-    if max < threshold:
-        diffs[:,ID] = 0
-        diffsum = 0
-        print('Detected differences of less than {} when comparing {} field component, therefore set as zero.'.format(threshold, fields[ID]))
-    else:
-        diffs[:,ID] = (np.abs(data1[:fieldssmallest,ID] - data2[:fieldssmallest,ID]) / max) * 100
-        diffsum = (np.sum(np.abs(data1[:fieldssmallest,ID] - data2[:fieldssmallest,ID])) / np.sum(np.abs(data1[:fieldssmallest,ID]))) * 100
-    print('Total differences in field component {}: {:.1f}%'.format(name, diffsum))
-
-# Plot data1
-fig1, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num=filename1 + ' versus ' + filename2, figsize=(20, 10), facecolor='w', edgecolor='w')
-ax1.plot(time1, data1[:,0],'r', lw=2, label='Ex')
-ax3.plot(time1, data1[:,2],'r', lw=2, label='Ey')
-ax5.plot(time1, data1[:,4],'r', lw=2, label='Ez')
-ax2.plot(time1, data1[:,1],'b', lw=2, label='Hx')
-ax4.plot(time1, data1[:,3],'b', lw=2, label='Hy')
-ax6.plot(time1, data1[:,5],'b', lw=2, label='Hz')
-    
-# Set ylabels
-ylabels = ['$E_x$, field strength [V/m]', '$H_x$, field strength [A/m]', '$E_y$, field strength [V/m]', '$H_y$, field strength [A/m]', '$E_z$, field strength [V/m]', '$H_z$, field strength [A/m]']
-[ax.set_ylabel(ylabels[index]) for index, ax in enumerate(fig1.axes)]
-
-# Turn on grid
-[ax.grid() for ax in fig1.axes]
-
-# Add data2 and set legend
-for index, ax in enumerate(fig1.axes):
-    if plotorder[index] in [0, 2, 4]:
-        ax.plot(time2, data2[:,plotorder[index]], 'r', label='data2', lw=2, ls='--')
-    else:
-        ax.plot(time2, data2[:,plotorder[index]], label='data2', lw=2, ls='--')
-    ax.set_xlim(0, time2[-1])
-    handles, existlabels = ax.get_legend_handles_labels()
-    ax.legend(handles, [os.path.split(filename1)[1], os.path.split(filename2)[1]])
-
-# Plots of differences
-fig2, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num='Deltas: ' + filename1 + ' versus ' + filename2, figsize=(20, 10), facecolor='w', edgecolor='w')
-ax1.plot(time1[:timesmallest], diffs[:,0],'r', lw=2, label='Ex')
-ax3.plot(time1[:timesmallest], diffs[:,2],'r', lw=2, label='Ey')
-ax5.plot(time1[:timesmallest], diffs[:,4],'r', lw=2, label='Ez')
-ax2.plot(time1[:timesmallest], diffs[:,1],'b', lw=2, label='Hx')
-ax4.plot(time1[:timesmallest], diffs[:,3],'b', lw=2, label='Hy')
-ax6.plot(time1[:timesmallest], diffs[:,5],'b', lw=2, label='Hz')
-
-# Set ylabels
-ylabels = ['$E_x$', '$H_x$', '$E_y$', '$H_y$', '$E_z$', '$H_z$']
-ylabels = [ylabel + ', percentage difference [%]' for ylabel in ylabels]
-[ax.set_ylabel(ylabels[index]) for index, ax in enumerate(fig2.axes)]
-
-# Set axes limits and turn on grid
-[ax.grid() for ax in fig2.axes]
-[ax.set_xlim(0, time1[timesmallest - 1]) for ax in fig2.axes]
-[ax.set_ylim(0, np.ceil(np.amax(np.abs(diffs)))) for ax in fig2.axes]
-
-# Show/print plots
-savename = os.path.abspath(os.path.dirname(filename1)) + os.sep + os.path.splitext(os.path.split(filename1)[1])[0] + '_vs_' + os.path.splitext(os.path.split(filename2)[1])[0]
-#fig1.savefig(savename + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-#fig2.savefig(savename + '_diffs.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-plt.show()
-

tests/test_models_basic.py

@@ -0,0 +1,232 @@
+# Copyright (C) 2015-2017: The University of Edinburgh
+#                 Authors: Craig Warren and Antonis Giannopoulos
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
+
+import datetime
+import os
+import sys
+from time import perf_counter
+
+from colorama import init, Fore, Style
+init()
+import h5py
+import numpy as np
+np.seterr(invalid='raise')
+import matplotlib.pyplot as plt
+
+if sys.platform == 'linux':
+    plt.switch_backend('agg')
+
+from gprMax.gprMax import api
+from gprMax.exceptions import GeneralError
+from tests.analytical_solutions import hertzian_dipole_fs
+
+"""Compare field outputs
+    
+    Usage:
+        cd gprMax
+        python -m tests.test_basic_suite
+    
+"""
+
+basepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'models_basic')
+
+# List of available test models
+testmodels = ['hertzian_dipole_fs_analytical', '2D_ExHyHz', '2D_EyHxHz', '2D_EzHxHy', 'cylinder_Ascan_2D', 'hertzian_dipole_fs', 'hertzian_dipole_hs', 'hertzian_dipole_dispersive']
+
+# Select a specific model if desired
+#testmodels = [testmodels[0], testmodels[1], testmodels[2], testmodels[3], testmodels[4], testmodels[5]]
+#testmodels = [testmodels[5]]
+testresults = dict.fromkeys(testmodels)
+path = '/rxs/rx1/'
+
+starttime = perf_counter()
+
+for i, model in enumerate(testmodels):
+    
+    testresults[model] = {}
+    
+    # Run model
+    api(os.path.join(basepath, model + os.path.sep + model + '.in'))
+
+    # Special case for analytical comparison
+    if model == 'hertzian_dipole_fs_analytical':
+        # Get output for model file
+        filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'),'r')
+        testresults[model]['Test version'] = filetest.attrs['gprMax']
+
+        # Get available field output component names
+        outputstest = list(filetest[path].keys())
+        
+        # Arrays for storing time
+        floattype = filetest[path + outputstest[0]].dtype
+        timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattype)
+        timetest = np.arange(0, filetest.attrs['dt'] * filetest.attrs['Iterations'], filetest.attrs['dt']) / 1e-9
+        timeref = timetest
+
+        # Arrays for storing field data
+        datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)), dtype=floattype)
+        print(datatest.shape)
+        for ID, name in enumerate(outputstest):
+            datatest[:,ID] = filetest[path + str(name)][:]
+
+        # Tx/Rx position to feed to analytical solution
+        rxpos = filetest[path].attrs['Position']
+        txpos = filetest['/srcs/src1/'].attrs['Position']
+        rxposrelative = ((rxpos[0] - txpos[0]), (rxpos[1] - txpos[1]), (rxpos[2] - txpos[2]))        
+
+        # Analytical solution of a dipole in free space
+        dataref = hertzian_dipole_fs(filetest.attrs['Iterations'], filetest.attrs['dt'], filetest.attrs['dx, dy, dz'], rxposrelative)
+
+        filetest.close()
+
+        # Diffs
+        datadiffs = np.zeros(datatest.shape, dtype=floattype)
+        for i in range(len(outputstest)):
+            max = np.amax(np.abs(dataref[:,i]))
+            try:
+                datadiffs[:,i] = ((np.abs(dataref[:,i] - datatest[:,i])) / max) * 100
+            except FloatingPointError:
+                print('FloatingPointError')
+                datadiffs[:,i] = 0
+
+        # Register test passed
+        threshold = 2 # Percent
+        if np.amax(np.amax(datadiffs)) < 2:
+            testresults[model]['Pass'] = True
+        else:
+            testresults[model]['Pass'] = False
+        testresults[model]['Max diff'] = np.amax(np.amax(datadiffs))
+
+    else:
+        # Get output for model and reference files
+        fileref = h5py.File(os.path.join(basepath, model + os.path.sep + model + '_ref.out'),'r')
+        filetest = h5py.File(os.path.join(basepath, model + os.path.sep + model + '.out'),'r')
+        testresults[model]['Ref version'] = fileref.attrs['gprMax']
+        testresults[model]['Test version'] = filetest.attrs['gprMax']
+        
+        # Get available field output component names
+        outputsref = list(fileref[path].keys())
+        outputstest = list(filetest[path].keys())
+        if outputsref != outputstest:
+            raise GeneralError('Field output components do not match reference solution')
+                             
+        # Check that type of float used to store fields matches
+        if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
+            raise GeneralError('Type of floating point number does not match reference solution')
+        else:
+            floattype = fileref[path + outputsref[0]].dtype
+
+        # Array for storing time
+        timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattype)
+        timeref = np.arange(0, fileref.attrs['dt'] * fileref.attrs['Iterations'], fileref.attrs['dt']) / 1e-9
+        timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattype)
+        timetest = np.arange(0, filetest.attrs['dt'] * filetest.attrs['Iterations'], filetest.attrs['dt']) / 1e-9
+
+        # Get available field output component names
+        outputsref = list(fileref[path].keys())
+        outputstest = list(filetest[path].keys())
+        if outputsref != outputstest:
+            raise GeneralError('Field output components do not match reference solution')
+
+        # Arrays for storing field data
+        dataref = np.zeros((fileref.attrs['Iterations'], len(outputsref)), dtype=floattype)
+        datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)), dtype=floattype)
+        for ID, name in enumerate(outputsref):
+            dataref[:,ID] = fileref[path + str(name)][:]
+            datatest[:,ID] = filetest[path + str(name)][:]
+
+        fileref.close()
+        filetest.close()
+
+        # Diffs
+        datadiffs = np.zeros(datatest.shape, dtype=floattype)
+        for i in range(len(outputstest)):
+            max = np.nanmax(np.abs(dataref[:,i]))
+            try:
+                datadiffs[:,i] = ((np.abs(dataref[:,i] - datatest[:,i])) / max) * 100
+            except FloatingPointError:
+                print('FloatingPointError')
+                datadiffs[:,i] = 0
+
+        # Register test passed
+        if not np.any(datadiffs):
+            testresults[model]['Pass'] = True
+        else:
+            testresults[model]['Pass'] = False
+        testresults[model]['Max diff'] = np.amax(np.amax(datadiffs))
+
+    # Plot datasets
+    fig1, ((ex1, hx1), (ey1, hy1), (ez1, hz1)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num=model + '.in', figsize=(20, 10), facecolor='w', edgecolor='w')
+    ex1.plot(timetest, datatest[:,0],'r', lw=2, label=model)
+    ex1.plot(timeref, dataref[:,0],'g', lw=2, ls='--', label=model + '(Ref)')
+    ey1.plot(timetest, datatest[:,1],'r', lw=2, label=model)
+    ey1.plot(timeref, dataref[:,1],'g', lw=2, ls='--', label=model + '(Ref)')
+    ez1.plot(timetest, datatest[:,2],'r', lw=2, label=model)
+    ez1.plot(timeref, dataref[:,2],'g', lw=2, ls='--', label=model + '(Ref)')
+    hx1.plot(timetest, datatest[:,3],'r', lw=2, label=model)
+    hx1.plot(timeref, dataref[:,3],'g', lw=2, ls='--', label=model + '(Ref)')
+    hy1.plot(timetest, datatest[:,4],'r', lw=2, label=model)
+    hy1.plot(timeref, dataref[:,4],'g', lw=2, ls='--', label=model + '(Ref)')
+    hz1.plot(timetest, datatest[:,5],'r', lw=2, label=model)
+    hz1.plot(timeref, dataref[:,5],'g', lw=2, ls='--', label=model + '(Ref)')
+    ylabels = ['$E_x$, field strength [V/m]', '$H_x$, field strength [A/m]', '$E_y$, field strength [V/m]', '$H_y$, field strength [A/m]', '$E_z$, field strength [V/m]', '$H_z$, field strength [A/m]']
+    for i, ax in enumerate(fig1.axes):
+        ax.set_ylabel(ylabels[i])
+        ax.set_xlim(0, np.amax(timetest))
+        ax.grid()
+        ax.legend()
+
+    # Plot diffs
+    fig2, ((ex2, hx2), (ey2, hy2), (ez2, hz2)) = plt.subplots(nrows=3, ncols=2, sharex=False, sharey='col', subplot_kw=dict(xlabel='Time [ns]'), num='Diffs: ' + model + '.in', figsize=(20, 10), facecolor='w', edgecolor='w')
+    ex2.plot(timeref, datadiffs[:,0],'r', lw=2, label='Ex')
+    ey2.plot(timeref, datadiffs[:,1],'r', lw=2, label='Ey')
+    ez2.plot(timeref, datadiffs[:,2],'r', lw=2, label='Ez')
+    hx2.plot(timeref, datadiffs[:,3],'r', lw=2, label='Hx')
+    hy2.plot(timeref, datadiffs[:,4],'r', lw=2, label='Hy')
+    hz2.plot(timeref, datadiffs[:,5],'r', lw=2, label='Hz')
+    ylabels = ['$E_x$, difference [%]', '$H_x$, difference [%]', '$E_y$, difference [%]', '$H_y$, difference [%]', '$E_z$, difference [%]', '$H_z$, difference [%]']
+    for i, ax in enumerate(fig2.axes):
+        ax.set_ylabel(ylabels[i])
+        ax.set_xlim(0, np.amax(timetest))
+        ax.grid()
+
+    # Save a PDF/PNG of the figure
+    savename = os.path.join(basepath, model + os.path.sep + model)
+    #fig1.savefig(savename + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+    #fig2.savefig(savename + '_diffs.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+    fig1.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
+    fig2.savefig(savename + '_diffs.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
+
+stoptime = perf_counter()
+
+# Summary of results
+passed = 0
+for name, data in testresults.items():
+    if 'analytical' in name:
+        if data['Pass']:
+            print(Fore.GREEN + "Test '{}.in' using v.{} compared to analytical solution passed. Maximum difference = {}%".format(name, data['Test version'], data['Max diff']) + Style.RESET_ALL)
+            passed += 1
+        else:
+            print(Fore.RED + "Test '{}.in' using v.{} compared to analytical solution failed. Maximum difference = {}%".format(name, data['Test version'], data['Max diff']) + Style.RESET_ALL)
+    else:
+        if data['Pass']:
+            print(Fore.GREEN + "Test '{}.in' using v.{} compared to reference solution using v.{} passed. Maximum difference = {}%".format(name, data['Test version'], data['Ref version'], data['Max diff']) + Style.RESET_ALL)
+            passed += 1
+        else:
+            print(Fore.RED + "Test '{}.in' using v.{} compared to reference solution using v.{} failed. Maximum difference = {}%".format(name, data['Test version'], data['Ref version'], data['Max diff']) + Style.RESET_ALL)
+print('{} of {} tests passed successfully in [HH:MM:SS]: {}'.format(passed, len(testmodels), datetime.timedelta(seconds=int(stoptime - starttime))))