diff --git a/testing/diff_output_files.py b/testing/diff_output_files.py
new file mode 100644
index 00000000..b2a3da6b
--- /dev/null
+++ b/testing/diff_output_files.py
@@ -0,0 +1,93 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# 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 logging
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+logger = logging.getLogger(__name__)
+
+def diff_output_files(filename1, filename2):
+    """Calculates differences between two output files.
+
+    Args:
+        filename1: string of filename (including path) of output file 1.
+        filename2: string of filename (including path) of output file 2.
+
+    Returns:
+        time: numpy array containing time.
+        datadiffs: numpy array containing power (dB) of differences.
+    """
+
+    file1 = h5py.File(Path(filename1), "r")
+    file2 = h5py.File(Path(filename2), "r")
+    # Path to receivers in files
+    path = "rxs/rx1/"
+
+    # Get available field output component names
+    outputs1 = list(file1[path].keys())
+    outputs2 = list(file2[path].keys())
+    if outputs1 != outputs2:
+        logger.exception("Field output components are not the same in each file")
+        raise ValueError
+
+    # Check that type of float used to store fields matches
+    floattype1 = file1[path + outputs1[0]].dtype
+    floattype2 = file2[path + outputs2[0]].dtype
+    if floattype1 != floattype2:
+        logger.warning(
+            f"Type of floating point number in test model ({file1[path + outputs1[0]].dtype}) "
+            f"does not match type in reference solution ({file2[path + outputs2[0]].dtype})\n"
+        )
+    
+    # Arrays for storing time
+    time1 = np.zeros((file1.attrs["Iterations"]), dtype=floattype1)
+    time1 = np.linspace(0, (file1.attrs["Iterations"] - 1), num=file1.attrs["Iterations"])
+    time2 = np.zeros((file2.attrs["Iterations"]), dtype=floattype2)
+    time2 = np.linspace(0, (file2.attrs["Iterations"] - 1), num=file2.attrs["Iterations"])
+
+    # Arrays for storing field data
+    data1 = np.zeros((file1.attrs["Iterations"], len(outputs1)), dtype=floattype1)
+    data2 = np.zeros((file2.attrs["Iterations"], len(outputs2)), dtype=floattype2)
+    for ID, name in enumerate(outputs1):
+        data1[:, ID] = file1[path + str(name)][:]
+        data2[:, ID] = file2[path + str(name)][:]
+        if np.any(np.isnan(data1[:, ID])) or np.any(np.isnan(data2[:, ID])):
+            logger.exception("Data contains NaNs")
+            raise ValueError
+
+    file1.close()
+    file2.close()
+
+    # Diffs
+    datadiffs = np.zeros(data1.shape, dtype=np.float64)
+    for i in range(len(outputs2)):
+        maxi = np.amax(np.abs(data1[:, i]))
+        datadiffs[:, i] = np.divide(
+            np.abs(data2[:, i] - data1[:, i]), maxi, out=np.zeros_like(data1[:, i]), where=maxi != 0
+        )  # Replace any division by zero with zero
+
+        # Calculate power (ignore warning from taking a log of any zero values)
+        with np.errstate(divide="ignore"):
+            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])
+        # Replace any NaNs or Infs from zero division
+        datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0
+
+    return time1, datadiffs
diff --git a/testing/models_pmls/plot_pml_comparison.py b/testing/models_pmls/plot_pml_comparison.py
deleted file mode 100644
index 95c5a79f..00000000
--- a/testing/models_pmls/plot_pml_comparison.py
+++ /dev/null
@@ -1,135 +0,0 @@
-# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
-#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
-#
-# 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 itertools
-import logging
-from operator import add
-from pathlib import Path
-
-import h5py
-import matplotlib.pyplot as plt
-import numpy as np
-
-logger = logging.getLogger(__name__)
-
-# Create/setup plot figure
-# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
-# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"]
-colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"]
-# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"]
-colors = itertools.cycle(colorIDs)
-# for i in range(2):
-#     next(colors)
-lines = itertools.cycle(("--", ":", "-.", "-"))
-markers = ["o", "d", "^", "s", "*"]
-
-parts = Path(__file__).parts
-path = "rxs/rx1/"
-basename = "pml_3D_pec_plate"
-PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"]
-maxerrors = []
-testmodels = ["pml_3D_pec_plate_" + s for s in PMLIDs]
-
-fig, ax = plt.subplots(
-    subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w"
-)
-
-for x, model in enumerate(testmodels):
-    # Open output file and read iterations
-    fileref = h5py.File(Path(*parts[:-1], basename, basename + "_ref.h5"), "r")
-    filetest = h5py.File(Path(*parts[:-1], basename, basename + str(x + 1) + ".h5"), "r")
-
-    # Get available field output component names
-    outputsref = list(fileref[path].keys())
-    outputstest = list(filetest[path].keys())
-    if outputsref != outputstest:
-        logger.exception("Field output components do not match reference solution")
-        raise ValueError
-
-    # Check that type of float used to store fields matches
-    if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
-        logger.warning(
-            f"Type of floating point number in test model ({filetest[path + outputstest[0]].dtype}) "
-            f"does not match type in reference solution ({fileref[path + outputsref[0]].dtype})\n"
-        )
-    floattyperef = fileref[path + outputsref[0]].dtype
-    floattypetest = filetest[path + outputstest[0]].dtype
-    # logger.info(f'Data type: {floattypetest}')
-
-    # Arrays for storing time
-    # timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattyperef)
-    # timeref = np.linspace(0, (fileref.attrs['Iterations'] - 1) * fileref.attrs['dt'], num=fileref.attrs['Iterations']) / 1e-9
-    # timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattypetest)
-    # timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'], num=filetest.attrs['Iterations']) / 1e-9
-    timeref = np.zeros((fileref.attrs["Iterations"]), dtype=floattyperef)
-    timeref = np.linspace(0, (fileref.attrs["Iterations"] - 1), num=fileref.attrs["Iterations"])
-    timetest = np.zeros((filetest.attrs["Iterations"]), dtype=floattypetest)
-    timetest = np.linspace(0, (filetest.attrs["Iterations"] - 1), num=filetest.attrs["Iterations"])
-
-    # Arrays for storing field data
-    dataref = np.zeros((fileref.attrs["Iterations"], len(outputsref)), dtype=floattyperef)
-    datatest = np.zeros((filetest.attrs["Iterations"], len(outputstest)), dtype=floattypetest)
-    for ID, name in enumerate(outputsref):
-        dataref[:, ID] = fileref[path + str(name)][:]
-        datatest[:, ID] = filetest[path + str(name)][:]
-        if np.any(np.isnan(datatest[:, ID])):
-            logger.exception("Test data contains NaNs")
-            raise ValueError
-
-    fileref.close()
-    filetest.close()
-
-    # Diffs
-    datadiffs = np.zeros(datatest.shape, dtype=np.float64)
-    for i in range(len(outputstest)):
-        maxi = np.amax(np.abs(dataref[:, i]))
-        datadiffs[:, i] = np.divide(
-            np.abs(datatest[:, i] - dataref[:, i]), maxi, out=np.zeros_like(dataref[:, i]), where=maxi != 0
-        )  # Replace any division by zero with zero
-
-        # Calculate power (ignore warning from taking a log of any zero values)
-        with np.errstate(divide="ignore"):
-            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])
-        # Replace any NaNs or Infs from zero division
-        datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0
-
-    # Print maximum error value
-    start = 210
-    maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]")
-    logger.info(f"{model}: Max. error {maxerrors[x]}")
-
-    # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
-    ax.plot(timeref[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model)
-    ax.set_xticks(np.arange(0, 2200, step=100))
-    ax.set_xlim([0, 2100])
-    ax.set_yticks(np.arange(-160, 0, step=20))
-    ax.set_ylim([-160, -20])
-    ax.set_axisbelow(True)
-    ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
-
-mylegend = list(map(add, PMLIDs, maxerrors))
-legend = ax.legend(mylegend, loc=1, fontsize=14)
-frame = legend.get_frame()
-frame.set_edgecolor("white")
-frame.set_alpha(0)
-
-plt.show()
-
-# Save a PDF/PNG of the figure
-# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
-# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
diff --git a/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py
new file mode 100644
index 00000000..52968ef1
--- /dev/null
+++ b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py
@@ -0,0 +1,80 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# 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 itertools
+import logging
+from operator import add
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from testing.diff_output_files import diff_output_files
+
+logger = logging.getLogger(__name__)
+
+# Create/setup plot figure
+# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
+# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"]
+colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"]
+# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"]
+colors = itertools.cycle(colorIDs)
+# for i in range(2):
+#     next(colors)
+lines = itertools.cycle(("--", ":", "-.", "-"))
+markers = ["o", "d", "^", "s", "*"]
+
+fn = Path(__file__)
+basename = "pml_3D_pec_plate"
+PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"]
+maxerrors = []
+testmodels = [basename + "_" + s for s in PMLIDs]
+
+fig, ax = plt.subplots(
+    subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w"
+)
+
+for x, model in enumerate(testmodels):
+    time, datadiffs = diff_output_files(fn.parent.joinpath(basename + "_ref.h5"), 
+                                        fn.parent.joinpath(basename + str(x + 1) + ".h5"))
+
+    # Print maximum error value
+    start = 210
+    maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]")
+    logger.info(f"{model}: Max. error {maxerrors[x]}")
+
+    # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
+    ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model)
+    ax.set_xticks(np.arange(0, 2200, step=100))
+    ax.set_xlim([0, 2100])
+    ax.set_yticks(np.arange(-160, 0, step=20))
+    ax.set_ylim([-160, -20])
+    ax.set_axisbelow(True)
+    ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
+
+mylegend = list(map(add, PMLIDs, maxerrors))
+legend = ax.legend(mylegend, loc=1, fontsize=14)
+frame = legend.get_frame()
+frame.set_edgecolor("white")
+frame.set_alpha(0)
+
+plt.show()
+
+# Save a PDF/PNG of the figure
+# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
diff --git a/testing/models_pmls/pml_basic/plot_pml_comparison.py b/testing/models_pmls/pml_basic/plot_pml_comparison.py
new file mode 100644
index 00000000..7429dca2
--- /dev/null
+++ b/testing/models_pmls/pml_basic/plot_pml_comparison.py
@@ -0,0 +1,75 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+#                 Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# 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 itertools
+import logging
+from operator import add
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from testing.diff_output_files import diff_output_files
+
+logger = logging.getLogger(__name__)
+
+# Create/setup plot figure
+# Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php
+colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"]
+colors = itertools.cycle(colorIDs)
+lines = itertools.cycle(("--", ":", "-.", "-"))
+markers = ["o", "d", "^", "s", "*"]
+
+fn = Path(__file__)
+basename = "pml_basic"
+PMLIDs = ["off", "x0", "y0", "z0", "xmax", "ymax", "zmax"]
+maxerrors = []
+
+for x in range(len(PMLIDs)):
+    file1 = fn.parent.joinpath(basename + str(x + 1) + "_CPU.h5")
+    file2 = fn.parent.joinpath(basename + str(x + 1) + "_GPU.h5")
+    time, datadiffs = diff_output_files(file1, file2)
+
+    # Print maximum error value
+    start = 0
+    maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]")
+    logger.info(f"{file1.name} - {file2.name}: Max. error {maxerrors[x]}")
+
+    fig, ax = plt.subplots(
+    subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w")
+
+    # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
+    ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=f"{file1.name} - {file2.name}")
+    ax.set_xticks(np.arange(0, 2200, step=100))
+    ax.set_xlim([0, 2100])
+    ax.set_yticks(np.arange(-160, 0, step=20))
+    ax.set_ylim([-160, -20])
+    ax.set_axisbelow(True)
+    ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed")
+
+mylegend = list(map(add, PMLIDs, maxerrors))
+legend = ax.legend(mylegend, loc=1, fontsize=14)
+frame = legend.get_frame()
+frame.set_edgecolor("white")
+frame.set_alpha(0)
+
+plt.show()
+
+# Save a PDF/PNG of the figure
+# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
+# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
diff --git a/testing/models_pmls/pml_basic.py b/testing/models_pmls/pml_basic/pml_basic.py
similarity index 100%
rename from testing/models_pmls/pml_basic.py
rename to testing/models_pmls/pml_basic/pml_basic.py