Renamed some modules.

gprMax/fields_updates.pyx

@@ -64,33 +64,33 @@ cpdef void update_electric(int nx, int ny, int nz, int nthreads, floattype_t[:,
 
     # 3D
     else:
-        for i in prange(0, nx - 1, nogil=True, schedule='static', num_threads=nthreads):
-            for j in range(0, ny - 1):
-                for k in range(0, nz - 1):
-                    materialEx = ID[0, i, j + 1, k + 1]
-                    materialEy = ID[1, i + 1, j, k + 1]
-                    materialEz = ID[2, i + 1, j + 1, k]
-                    Ex[i, j + 1, k + 1] = updatecoeffsE[materialEx, 0] * Ex[i, j + 1, k + 1] + updatecoeffsE[materialEx, 2] * (Hz[i, j + 1, k + 1] - Hz[i, j, k + 1]) - updatecoeffsE[materialEx, 3] * (Hy[i, j + 1, k + 1] - Hy[i, j + 1, k])
-                    Ey[i + 1, j, k + 1] = updatecoeffsE[materialEy, 0] * Ey[i + 1, j, k + 1] + updatecoeffsE[materialEy, 3] * (Hx[i + 1, j, k + 1] - Hx[i + 1, j, k]) - updatecoeffsE[materialEy, 1] * (Hz[i + 1, j, k + 1] - Hz[i, j, k + 1])
-                    Ez[i + 1, j + 1, k] = updatecoeffsE[materialEz, 0] * Ez[i + 1, j + 1, k] + updatecoeffsE[materialEz, 1] * (Hy[i + 1, j + 1, k] - Hy[i, j + 1, k]) - updatecoeffsE[materialEz, 2] * (Hx[i + 1, j + 1, k] - Hx[i + 1, j, k])
-
-        # Ex components at nx - 1
-        for j in prange(1, ny, nogil=True, schedule='static', num_threads=nthreads):
-            for k in range(1, nz):
-                materialEx = ID[0, nx - 1, j, k]
-                Ex[nx - 1, j, k] = updatecoeffsE[materialEx, 0] * Ex[nx - 1, j, k] + updatecoeffsE[materialEx, 2] * (Hz[nx - 1, j, k] - Hz[nx - 1, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[nx - 1, j, k] - Hy[nx - 1, j, k - 1])
-
-        # Ey components at ny - 1
-        for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
-            for k in range(1, nz):
-                materialEy = ID[1, i, ny - 1, k]
-                Ey[i, ny - 1, k] = updatecoeffsE[materialEy, 0] * Ey[i, ny - 1, k] + updatecoeffsE[materialEy, 3] * (Hx[i, ny - 1, k] - Hx[i, ny - 1, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, ny - 1, k] - Hz[i - 1, ny - 1, k])
-
-        # Ez components at nz - 1
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
-                materialEz = ID[2, i, j, nz - 1]
-                Ez[i, j, nz - 1] = updatecoeffsE[materialEz, 0] * Ez[i, j, nz - 1] + updatecoeffsE[materialEz, 1] * (Hy[i, j, nz - 1] - Hy[i - 1, j, nz - 1]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, nz - 1] - Hx[i, j - 1, nz - 1])
+                for k in range(1, nz):
+                    materialEx = ID[0, i, j, k]
+                    materialEy = ID[1, i, j, k]
+                    materialEz = ID[2, i, j, k]
+                    Ex[i, j, k] = updatecoeffsE[materialEx, 0] * Ex[i, j, k] + updatecoeffsE[materialEx, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[i, j, k] - Hy[i, j, k - 1])
+                    Ey[i, j, k] = updatecoeffsE[materialEy, 0] * Ey[i, j, k] + updatecoeffsE[materialEy, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, j, k] - Hz[i - 1, j, k])
+                    Ez[i, j, k] = updatecoeffsE[materialEz, 0] * Ez[i, j, k] + updatecoeffsE[materialEz, 1] * (Hy[i, j, k] - Hy[i - 1, j, k]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, k] - Hx[i, j - 1, k])
+
+        # Ex components at i = 0
+        for j in prange(1, ny, nogil=True, schedule='static', num_threads=nthreads):
+            for k in range(1, nz):
+                materialEx = ID[0, 0, j, k]
+                Ex[0, j, k] = updatecoeffsE[materialEx, 0] * Ex[0, j, k] + updatecoeffsE[materialEx, 2] * (Hz[0, j, k] - Hz[0, j - 1, k]) - updatecoeffsE[materialEx, 3] * (Hy[0, j, k] - Hy[0, j, k - 1])
+
+        # Ey components at j = 0
+        for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
+            for k in range(1, nz):
+                materialEy = ID[1, i, 0, k]
+                Ey[i, 0, k] = updatecoeffsE[materialEy, 0] * Ey[i, 0, k] + updatecoeffsE[materialEy, 3] * (Hx[i, 0, k] - Hx[i, 0, k - 1]) - updatecoeffsE[materialEy, 1] * (Hz[i, 0, k] - Hz[i - 1, 0, k])
+
+        # Ez components at k = 0
+        for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
+            for j in range(1, ny):
+                materialEz = ID[2, i, j, 0]
+                Ez[i, j, 0] = updatecoeffsE[materialEz, 0] * Ez[i, j, 0] + updatecoeffsE[materialEz, 1] * (Hy[i, j, 0] - Hy[i - 1, j, 0]) - updatecoeffsE[materialEz, 2] * (Hx[i, j, 0] - Hx[i, j - 1, 0])
 
 
 #################################################
@@ -111,9 +111,7 @@ cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nt
     cdef float phi = 0
 
     # Ex component
-    if ny == 1 or nz == 1:
-        pass
-    else:
+    if ny != 1 or nz != 1:
         for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(1, nz):
@@ -125,9 +123,7 @@ cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nt
                     Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
 
     # Ey component
-    if nx == 1 or nz == 1:
-        pass
-    else:
+    if nx != 1 or nz != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(0, ny):
                 for k in range(1, nz):
@@ -139,9 +135,7 @@ cpdef void update_electric_dispersive_multipole_A(int nx, int ny, int nz, int nt
                     Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
 
     # Ez component
-    if nx == 1 or ny == 1:
-        pass
-    else:
+    if nx != 1 or ny != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(0, nz):
@@ -168,9 +162,7 @@ cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nt
     cdef int material
 
     # Ex component
-    if ny == 1 or nz == 1:
-        pass
-    else:
+    if ny != 1 or nz != 1:
         for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(1, nz):
@@ -179,9 +171,7 @@ cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nt
                         Tx[pole, i, j, k] = Tx[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ex[i, j, k]
 
     # Ey component
-    if nx == 1 or nz == 1:
-        pass
-    else:
+    if nx != 1 or nz != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(0, ny):
                 for k in range(1, nz):
@@ -190,9 +180,7 @@ cpdef void update_electric_dispersive_multipole_B(int nx, int ny, int nz, int nt
                         Ty[pole, i, j, k] = Ty[pole, i, j, k] - updatecoeffsdispersive[material, 2 + (pole * 3)] * Ey[i, j, k]
 
     # Ez component
-    if nx == 1 or ny == 1:
-        pass
-    else:
+    if nx != 1 or ny != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(0, nz):
@@ -215,9 +203,7 @@ cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthrea
     cdef float phi = 0
 
     # Ex component
-    if ny == 1 or nz == 1:
-        pass
-    else:
+    if ny != 1 or nz != 1:
         for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(1, nz):
@@ -227,9 +213,7 @@ cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthrea
                     Ex[i, j, k] = updatecoeffsE[material, 0] * Ex[i, j, k] + updatecoeffsE[material, 2] * (Hz[i, j, k] - Hz[i, j - 1, k]) - updatecoeffsE[material, 3] * (Hy[i, j, k] - Hy[i, j, k - 1]) - updatecoeffsE[material, 4] * phi
 
     # Ey component
-    if nx == 1 or nz == 1:
-        pass
-    else:
+    if nx != 1 or nz != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(0, ny):
                 for k in range(1, nz):
@@ -239,9 +223,7 @@ cpdef void update_electric_dispersive_1pole_A(int nx, int ny, int nz, int nthrea
                     Ey[i, j, k] = updatecoeffsE[material, 0] * Ey[i, j, k] + updatecoeffsE[material, 3] * (Hx[i, j, k] - Hx[i, j, k - 1]) - updatecoeffsE[material, 1] * (Hz[i, j, k] - Hz[i - 1, j, k]) - updatecoeffsE[material, 4] * phi
 
     # Ez component
-    if nx == 1 or ny == 1:
-        pass
-    else:
+    if nx != 1 or ny != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(0, nz):
@@ -264,9 +246,7 @@ cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthrea
     cdef int material
 
     # Ex component
-    if ny == 1 or nz == 1:
-        pass
-    else:
+    if ny != 1 or nz != 1:
         for i in prange(0, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(1, nz):
@@ -274,9 +254,7 @@ cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthrea
                     Tx[0, i, j, k] = Tx[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ex[i, j, k]
 
     # Ey component
-    if nx == 1 or nz == 1:
-        pass
-    else:
+    if nx != 1 or nz != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(0, ny):
                 for k in range(1, nz):
@@ -284,9 +262,7 @@ cpdef void update_electric_dispersive_1pole_B(int nx, int ny, int nz, int nthrea
                     Ty[0, i, j, k] = Ty[0, i, j, k] - updatecoeffsdispersive[material, 2] * Ey[i, j, k]
 
     # Ez component
-    if nx == 1 or ny == 1:
-        pass
-    else:
+    if nx != 1 or ny != 1:
         for i in prange(1, nx, nogil=True, schedule='static', num_threads=nthreads):
             for j in range(1, ny):
                 for k in range(0, nz):

gprMax/model_build_run.py

@@ -0,0 +1,323 @@
+# 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 itertools
+import os
+import psutil
+import sys
+from time import perf_counter
+
+from colorama import init, Fore, Style
+init()
+import numpy as np
+from terminaltables import AsciiTable
+from tqdm import tqdm
+
+from gprMax.constants import floattype, cfloattype, ccomplextype
+from gprMax.exceptions import GeneralError
+from gprMax.fields_outputs import store_outputs, write_hdf5_outputfile
+from gprMax.fields_updates import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B
+from gprMax.grid import FDTDGrid, dispersion_analysis
+from gprMax.input_cmds_geometry import process_geometrycmds
+from gprMax.input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names
+from gprMax.input_cmds_multiuse import process_multicmds
+from gprMax.input_cmds_singleuse import process_singlecmds
+from gprMax.materials import Material, process_materials
+from gprMax.pml import PML, build_pmls
+from gprMax.utilities import get_terminal_width, human_size, open_path_file
+from gprMax.yee_cell_build import build_electric_components, build_magnetic_components
+
+
+def run_model(args, currentmodelrun, numbermodelruns, inputfile, usernamespace):
+    """Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop.
+
+    Args:
+        args (dict): Namespace with command line arguments
+        currentmodelrun (int): Current model run number.
+        numbermodelruns (int): Total number of model runs.
+        inputfile (str/object): Name or file object of the input file.
+        usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
+
+    Returns:
+        tsolve (int): Length of time (seconds) of main FDTD calculations
+    """
+
+    # Monitor memory usage
+    p = psutil.Process()
+
+    # Declare variable to hold FDTDGrid class
+    global G
+
+    # Normal model reading/building process; bypassed if geometry information to be reused
+    if 'G' not in globals():
+
+        # Initialise an instance of the FDTDGrid class
+        G = FDTDGrid()
+
+        with open_path_file(inputfile) as f:
+            G.inputfilename = os.path.split(f.name)[1]
+            G.inputdirectory = os.path.dirname(os.path.abspath(f.name))
+            inputfilestr = '\n--- Model {}/{}, input file: {}'.format(currentmodelrun, numbermodelruns, f.name)
+            print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
+
+            # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
+            usernamespace['current_model_run'] = currentmodelrun
+
+            # Read input file and process any Python or include commands
+            processedlines = process_python_include_code(f, usernamespace)
+
+        # Print constants/variables in user-accessable namespace
+        uservars = ''
+        for key, value in sorted(usernamespace.items()):
+            if key != '__builtins__':
+                uservars += '{}: {}, '.format(key, value)
+        print('Constants/variables used/available for Python scripting: {{{}}}\n'.format(uservars[:-2]))
+
+        # Write a file containing the input commands after Python or include commands have been processed
+        if args.write_processed:
+            write_processed_file(os.path.join(G.inputdirectory, G.inputfilename), currentmodelrun, numbermodelruns, processedlines)
+
+        # Check validity of command names and that essential commands are present
+        singlecmds, multicmds, geometry = check_cmd_names(processedlines)
+
+        # Create built-in materials
+        m = Material(0, 'pec')
+        m.se = float('inf')
+        m.type = 'builtin'
+        m.averagable = False
+        G.materials.append(m)
+        m = Material(1, 'free_space')
+        m.type = 'builtin'
+        G.materials.append(m)
+
+        # Process parameters for commands that can only occur once in the model
+        process_singlecmds(singlecmds, G)
+
+        # Process parameters for commands that can occur multiple times in the model
+        print()
+        process_multicmds(multicmds, G)
+
+        # Initialise an array for volumetric material IDs (solid), boolean arrays for specifying materials not to be averaged (rigid),
+        # an array for cell edge IDs (ID)
+        G.initialise_geometry_arrays()
+
+        # Initialise arrays for the field components
+        G.initialise_field_arrays()
+
+        # Process geometry commands in the order they were given
+        process_geometrycmds(geometry, G)
+
+        # Build the PMLs and calculate initial coefficients
+        print()
+        if all(value == 0 for value in G.pmlthickness.values()):
+            if G.messages:
+                print('PML boundaries: switched off')
+            pass  # If all the PMLs are switched off don't need to build anything
+        else:
+            if G.messages:
+                if all(value == G.pmlthickness['x0'] for value in G.pmlthickness.values()):
+                    pmlinfo = str(G.pmlthickness['x0']) + ' cells'
+                else:
+                    pmlinfo = ''
+                    for key, value in G.pmlthickness.items():
+                        pmlinfo += '{}: {} cells, '.format(key, value)
+                    pmlinfo = pmlinfo[:-2]
+                print('PML boundaries: {}'.format(pmlinfo))
+            pbar = tqdm(total=sum(1 for value in G.pmlthickness.values() if value > 0), desc='Building PML boundaries', ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
+            build_pmls(G, pbar)
+            pbar.close()
+
+        # Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
+        print()
+        pbar = tqdm(total=2, desc='Building main grid', ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
+        build_electric_components(G.solid, G.rigidE, G.ID, G)
+        pbar.update()
+        build_magnetic_components(G.solid, G.rigidH, G.ID, G)
+        pbar.update()
+        pbar.close()
+
+        # Process any voltage sources (that have resistance) to create a new material at the source location
+        for voltagesource in G.voltagesources:
+            voltagesource.create_material(G)
+
+        # Initialise arrays of update coefficients to pass to update functions
+        G.initialise_std_update_coeff_arrays()
+
+        # Initialise arrays of update coefficients and temporary values if there are any dispersive materials
+        if Material.maxpoles != 0:
+            G.initialise_dispersive_arrays()
+
+        # Process complete list of materials - calculate update coefficients, store in arrays, and build text list of materials/properties
+        materialsdata = process_materials(G)
+        if G.messages:
+            materialstable = AsciiTable(materialsdata)
+            materialstable.outer_border = False
+            materialstable.justify_columns[0] = 'right'
+            print(materialstable.table)
+
+        # Check to see if numerical dispersion might be a problem
+        results = dispersion_analysis(G)
+        if results['N'] < G.mingridsampling:
+            raise GeneralError("Non-physical wave propagation: Material '{}' has wavelength sampled by {} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {:g}Hz".format(results['material'].ID, results['N'], results['maxfreq']))
+        elif results['deltavp'] and np.abs(results['deltavp']) > G.maxnumericaldisp and G.messages:
+            print(Fore.RED + "\nWARNING: Potentially significant numerical dispersion. Estimated largest physical phase-velocity error is {:.2f}% in material '{}' whose wavelength sampled by {} cells. Maximum significant frequency estimated as {:g}Hz".format(results['deltavp'], results['material'].ID, results['N'], results['maxfreq']) + Style.RESET_ALL)
+        elif results['deltavp'] and G.messages:
+            print("\nNumerical dispersion analysis: estimated largest physical phase-velocity error is {:.2f}% in material '{}' whose wavelength sampled by {} cells. Maximum significant frequency estimated as {:g}Hz".format(results['deltavp'], results['material'].ID, results['N'], results['maxfreq']))
+
+    # If geometry information to be reused between model runs
+    else:
+        with open_path_file(inputfile) as f:
+            inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(currentmodelrun, numbermodelruns, f.name)
+            print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
+
+        # Clear arrays for field components
+        G.initialise_field_arrays()
+
+        # Clear arrays for fields in PML
+        for pml in G.pmls:
+            pml.initialise_field_arrays()
+
+    # Adjust position of simple sources and receivers if required
+    if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0:
+        for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
+            if currentmodelrun == 1:
+                if source.xcoord + G.srcsteps[0] * (numbermodelruns - 1) < 0 or source.xcoord + G.srcsteps[0] * (numbermodelruns - 1) > G.nx or source.ycoord + G.srcsteps[1] * (numbermodelruns - 1) < 0 or source.ycoord + G.srcsteps[1] * (numbermodelruns - 1) > G.ny or source.zcoord + G.srcsteps[2] * (numbermodelruns - 1) < 0 or source.zcoord + G.srcsteps[2] * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Source(s) will be stepped to a position outside the domain.')
+            source.xcoord = source.xcoordorigin + (currentmodelrun - 1) * G.srcsteps[0]
+            source.ycoord = source.ycoordorigin + (currentmodelrun - 1) * G.srcsteps[1]
+            source.zcoord = source.zcoordorigin + (currentmodelrun - 1) * G.srcsteps[2]
+    if G.rxsteps[0] != 0 or G.rxsteps[1] != 0 or G.rxsteps[2] != 0:
+        for receiver in G.rxs:
+            if currentmodelrun == 1:
+                if receiver.xcoord + G.rxsteps[0] * (numbermodelruns - 1) < 0 or receiver.xcoord + G.rxsteps[0] * (numbermodelruns - 1) > G.nx or receiver.ycoord + G.rxsteps[1] * (numbermodelruns - 1) < 0 or receiver.ycoord + G.rxsteps[1] * (numbermodelruns - 1) > G.ny or receiver.zcoord + G.rxsteps[2] * (numbermodelruns - 1) < 0 or receiver.zcoord + G.rxsteps[2] * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
+            receiver.xcoord = receiver.xcoordorigin + (currentmodelrun - 1) * G.rxsteps[0]
+            receiver.ycoord = receiver.ycoordorigin + (currentmodelrun - 1) * G.rxsteps[1]
+            receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2]
+
+    # Write files for any geometry views and geometry object outputs
+    if not (G.geometryviews or G.geometryobjectswrite) and args.geometry_only:
+        print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL)
+    if G.geometryviews:
+        print()
+        for i, geometryview in enumerate(G.geometryviews):
+            geometryview.set_filename(currentmodelrun, numbermodelruns, G)
+            pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
+            geometryview.write_vtk(currentmodelrun, numbermodelruns, G, pbar)
+            pbar.close()
+    if G.geometryobjectswrite:
+        for i, geometryobject in enumerate(G.geometryobjectswrite):
+            pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry object file {}/{}, {}'.format(i + 1, len(G.geometryobjectswrite), os.path.split(geometryobject.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
+            geometryobject.write_hdf5(G, pbar)
+            pbar.close()
+
+    # Run simulation (if not only looking ar geometry information)
+    if not args.geometry_only:
+
+        # Prepare any snapshot files
+        for snapshot in G.snapshots:
+            snapshot.prepare_vtk_imagedata(currentmodelrun, numbermodelruns, G)
+
+        # Output filename
+        inputfileparts = os.path.splitext(os.path.join(G.inputdirectory, G.inputfilename))
+        if numbermodelruns == 1:
+            outputfile = inputfileparts[0] + '.out'
+        else:
+            outputfile = inputfileparts[0] + str(currentmodelrun) + '.out'
+        print('\nOutput file: {}\n'.format(outputfile))
+
+        # Main FDTD solving functions for either CPU or GPU
+        tsolve = solve_cpu(currentmodelrun, numbermodelruns, G)
+
+        # Write an output file in HDF5 format
+        write_hdf5_outputfile(outputfile, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
+
+        if G.messages:
+            print('Memory (RAM) used: ~{}'.format(human_size(p.memory_info().rss)))
+            print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve)))
+
+        # If geometry information to be reused between model runs then FDTDGrid class instance must be global so that it persists
+        if not args.geometry_fixed:
+            del G
+
+        return tsolve
+
+
+def solve_cpu(currentmodelrun, numbermodelruns, G):
+    """Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for electric and magnetic field updates, and PML updates.
+
+    Args:
+        currentmodelrun (int): Current model run number.
+        numbermodelruns (int): Total number of model runs.
+        G (class): Grid class instance - holds essential parameters describing the model.
+
+    Returns:
+        tsolve (float): Time taken to execute solving
+    """
+
+    tsolvestart = perf_counter()
+
+    for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(numbermodelruns), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable):
+        # Store field component values for every receiver and transmission line
+        store_outputs(iteration, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
+
+        # Write any snapshots to file
+        for i, snap in enumerate(G.snapshots):
+            if snap.time == iteration + 1:
+                snapiters = 36 * (((snap.xf - snap.xs) / snap.dx) * ((snap.yf - snap.ys) / snap.dy) * ((snap.zf - snap.zs) / snap.dz))
+                pbar = tqdm(total=snapiters, leave=False, unit='byte', unit_scale=True, desc='  Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
+                snap.write_vtk_imagedata(G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G, pbar)
+                pbar.close()
+    
+        # Update magnetic field components
+        update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+        # Update magnetic field components with the PML correction
+        for pml in G.pmls:
+            pml.update_magnetic(G)
+
+        # Update magnetic field components from sources
+        for source in G.transmissionlines + G.magneticdipoles:
+            source.update_magnetic(iteration, G.updatecoeffsH, G.ID, G.Hx, G.Hy, G.Hz, G)
+
+        # Update electric field components
+        if Material.maxpoles == 0:  # All materials are non-dispersive so do standard update
+            update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+        elif Material.maxpoles == 1:  # If there are any dispersive materials do 1st part of dispersive update (it is split into two parts as it requires present and updated electric field values).
+            update_electric_dispersive_1pole_A(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+        elif Material.maxpoles > 1:
+            update_electric_dispersive_multipole_A(G.nx, G.ny, G.nz, G.nthreads, Material.maxpoles, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+        # Update electric field components with the PML correction
+        for pml in G.pmls:
+            pml.update_electric(G)
+
+        # Update electric field components from sources (update any Hertzian dipole sources last)
+        for source in G.voltagesources + G.transmissionlines + G.hertziandipoles:
+            source.update_electric(iteration, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G)
+
+        # If there are any dispersive materials do 2nd part of dispersive update (it is split into two parts as it requires present and updated electric field values). Therefore it can only be completely updated after the electric field has been updated by the PML and source updates.
+        if Material.maxpoles == 1:
+            update_electric_dispersive_1pole_B(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
+        elif Material.maxpoles > 1:
+            update_electric_dispersive_multipole_B(G.nx, G.ny, G.nz, G.nthreads, Material.maxpoles, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
+
+    tsolve = perf_counter() - tsolvestart
+
+    return tsolve

tests/models_basic/magnetic_dipole_fs/magnetic_dipole_fs.in

@@ -0,0 +1,8 @@
+#title: Magnetic 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
+#magnetic_dipole: z 0.050 0.050 0.050 myWave
+#rx: 0.070 0.070 0.070

tests/test_models.py

@@ -16,16 +16,13 @@
 # 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(divide='raise')
 import matplotlib.pyplot as plt
 
 if sys.platform == 'linux':
@@ -39,16 +36,21 @@ from tests.analytical_solutions import hertzian_dipole_fs
 
     Usage:
         cd gprMax
-        python -m tests.test_models_basic
+        python -m tests.test_models
 """
 
-basepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'models_basic')
+basepath = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'models_')
+basepath += 'basic'
+#basepath += 'advanced'
 
-# 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']
+# List of available basic 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', 'magnetic_dipole_fs']
+
+# List of available advanced test models
+#testmodels = ['antenna_GSSI_1500_fs', 'antenna_MALA_1200_fs']
 
 # Select a specific model if desired
-#testmodels = [testmodels[0], testmodels[1], testmodels[2], testmodels[3], testmodels[4], testmodels[5]]
+#testmodels = [testmodels[0], testmodels[5], testmodels[7]]
 #testmodels = [testmodels[5]]
 testresults = dict.fromkeys(testmodels)
 path = '/rxs/rx1/'
@@ -56,8 +58,6 @@ path = '/rxs/rx1/'
 # Minimum value of difference to plot (dB)
 plotmin = -160
 
-starttime = perf_counter()
-
 for i, model in enumerate(testmodels):
 
     testresults[model] = {}
@@ -97,10 +97,6 @@ for i, model in enumerate(testmodels):
 
         filetest.close()
 
-        # Threshold below which test is considered passed (dB)
-        threshold = -35
-        testresults[model]['Threshold'] = threshold
-
     else:
         # Get output for model and reference files
         fileref = h5py.File(os.path.join(basepath, model + os.path.sep + model + '_ref.out'), 'r')
@@ -116,25 +112,19 @@ for i, model in enumerate(testmodels):
 
         # 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
+            print(Fore.RED + 'WARNING: Type of floating point number in test model ({}) does not match type in reference solution ({})\n'.format(filetest[path + outputstest[0]].dtype, fileref[path + outputsref[0]].dtype) + Style.RESET_ALL)
+        floattyperef = fileref[path + outputsref[0]].dtype
+        floattypetest = filetest[path + outputstest[0]].dtype
 
         # Array for storing time
-        timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattype)
+        timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattyperef)
         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.zeros((filetest.attrs['Iterations']), dtype=floattypetest)
         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)
+        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)][:]
@@ -144,26 +134,16 @@ for i, model in enumerate(testmodels):
         fileref.close()
         filetest.close()
 
-        # Threshold below which test is considered passed (dB)
-        threshold = -120
-        testresults[model]['Threshold'] = threshold
-
     # Diffs
-    datadiffs = np.zeros(datatest.shape, dtype=floattype)
+    datadiffs = np.zeros(datatest.shape, dtype=np.float64)
     for i in range(len(outputstest)):
         max = np.amax(np.abs(dataref[:, i]))
-        try:
-            datadiffs[:, i] = 20 * np.log10(((np.abs(dataref[:, i] - datatest[:, i])) / max))
-        # If a divide by zero error is encountered, consider the difference to be minimum plotted
-        except FloatingPointError:
-            datadiffs[:, i] = plotmin
+        datadiffs[:, i] = np.divide(np.abs(dataref[:, i] - datatest[:, i]), max, out=np.zeros_like(dataref[:, i]), where=max!=0) # Replace any division by zero with zero
+        with np.errstate(divide='ignore'):
+            datadiffs[:, i] = 20 * np.log10(datadiffs[:, i]) # Ignore any zero division in log10
 
-    # Register test passed/failed
+    # Store max difference
     maxdiff = np.amax(np.amax(datadiffs))
-    if maxdiff <= threshold:
-        testresults[model]['Pass'] = True
-    else:
-        testresults[model]['Pass'] = False
     testresults[model]['Max diff'] = maxdiff
 
     # Plot datasets
@@ -199,7 +179,7 @@ for i, model in enumerate(testmodels):
     for i, ax in enumerate(fig2.axes):
         ax.set_ylabel(ylabels[i])
         ax.set_xlim(0, np.amax(timetest))
-        ax.set_ylim([plotmin, 0])
+        ax.set_ylim([plotmin, np.amax(np.amax(datadiffs))])
         ax.grid()
 
     # Save a PDF/PNG of the figure
@@ -209,21 +189,9 @@ for i, model in enumerate(testmodels):
     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():
+for name, data in sorted(testresults.items()):
     if 'analytical' in name:
-        if data['Pass']:
-            print(Fore.GREEN + "Test '{}.in' using v.{} compared to analytical solution passed. Max difference {:.2f}dB <= {:.2f}dB threshold".format(name, data['Test version'], data['Max diff'], data['Threshold']) + Style.RESET_ALL)
-            passed += 1
-        else:
-            print(Fore.RED + "Test '{}.in' using v.{} compared to analytical solution failed. Max difference {:.2f}dB <= {:.2f}dB threshold".format(name, data['Test version'], data['Max diff'], data['Threshold']) + Style.RESET_ALL)
+        print(Fore.CYAN + "Test '{}.in' using v.{} compared to analytical solution. Max difference {:.2f}dB.".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. Max difference {:.2f}dB <= {:.2f}dB threshold".format(name, data['Test version'], data['Ref version'], data['Max diff'], data['Threshold']) + Style.RESET_ALL)
-            passed += 1
-        else:
-            print(Fore.RED + "Test '{}.in' using v.{} compared to reference solution using v.{} failed. Max difference {:.2f}dB <= {:.2f}dB threshold".format(name, data['Test version'], data['Ref version'], data['Max diff'], data['Threshold']) + Style.RESET_ALL)
-print('{} of {} tests passed successfully in [HH:MM:SS]: {}'.format(passed, len(testmodels), datetime.timedelta(seconds=int(stoptime - starttime))))
+        print(Fore.CYAN + "Test '{}.in' using v.{} compared to reference solution using v.{}. Max difference {:.2f}dB.".format(name, data['Test version'], data['Ref version'], data['Max diff']) + Style.RESET_ALL)