Merge branch 'mpi' into add_reframe_tests

gprMax/cmds_multiuse.py

@@ -871,7 +871,7 @@ class Rx(UserObjectMulti):
 
         try:
             r.ID = self.kwargs["id"]
-            outputs = [self.kwargs["outputs"]]
+            outputs = list(self.kwargs["outputs"])
         except KeyError:
             # If no ID or outputs are specified, use default
             r.ID = f"{r.__class__.__name__}({str(r.xcoord)},{str(r.ycoord)},{str(r.zcoord)})"

gprMax/config.py

@@ -20,11 +20,14 @@ import logging
 import sys
 import warnings
 from pathlib import Path
+from typing import List, Optional, Union
 
 import cython
 import numpy as np
 from colorama import Fore, Style, init
 
+from gprMax.scene import Scene
+
 init()
 from scipy.constants import c
 from scipy.constants import epsilon_0 as e0
@@ -35,34 +38,17 @@ from .utilities.utilities import get_terminal_width
 
 logger = logging.getLogger(__name__)
 
-# Single instance of SimConfig to hold simulation configuration parameters.
-sim_config = None
-
-# Instances of ModelConfig that hold model configuration parameters.
-model_configs = []
-
-# Each model in a simulation is given a unique number when the instance of
-# ModelConfig is created
-model_num = 0
-
-
-def get_model_config():
-    """Return ModelConfig instace for specific model."""
-    if sim_config.args.taskfarm:
-        return model_configs
-    else:
-        return model_configs[model_num]
-
 
 class ModelConfig:
     """Configuration parameters for a model.
     N.B. Multiple models can exist within a simulation
     """
 
-    def __init__(self):
+    def __init__(self, model_num):
         self.mode = "3D"
         self.grids = []
         self.ompthreads = None
+        self.model_num = model_num
 
         # Store information for CUDA or OpenCL solver
         #   dev: compute device object.
@@ -73,7 +59,7 @@ class ModelConfig:
         if sim_config.general["solver"] in ["cuda", "opencl"]:
             if sim_config.general["solver"] == "cuda":
                 devs = sim_config.args.gpu
-            elif sim_config.general["solver"] == "opencl":
+            else:  # opencl
                 devs = sim_config.args.opencl
 
             # If a list of lists of deviceIDs is found, flatten it
@@ -87,20 +73,29 @@ class ModelConfig:
             except:
                 deviceID = 0
 
-            self.device = {"dev": sim_config.set_model_device(deviceID), "deviceID": deviceID, "snapsgpu2cpu": False}
+            self.device = {
+                "dev": sim_config.get_model_device(deviceID),
+                "deviceID": deviceID,
+                "snapsgpu2cpu": False,
+            }
 
         # Total memory usage for all grids in the model. Starts with 50MB overhead.
         self.mem_overhead = 65e6
         self.mem_use = self.mem_overhead
 
-        self.reuse_geometry = False
-
         # String to print at start of each model run
-        s = f"\n--- Model {model_num + 1}/{sim_config.model_end}, " f"input file: {sim_config.input_file_path}"
+        s = (
-        self.inputfilestr = Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
+            f"\n--- Model {model_num + 1}/{sim_config.model_end}, "
+            f"input file: {sim_config.input_file_path}"
+        )
+        self.inputfilestr = (
+            Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
+        )
 
         # Output file path and name for specific model
-        self.appendmodelnumber = "" if sim_config.args.n == 1 else str(model_num + 1)  # Indexed from 1
+        self.appendmodelnumber = (
+            "" if sim_config.args.n == 1 else str(model_num + 1)
+        )  # Indexed from 1
         self.set_output_file_path()
 
         # Numerical dispersion analysis parameters
@@ -128,17 +123,17 @@ class ModelConfig:
             "crealfunc": None,
         }
 
+    def reuse_geometry(self):
+        return self.model_num != 0 and sim_config.args.geometry_fixed
+
     def get_scene(self):
-        try:
+        return sim_config.get_scene(self.model_num)
-            return sim_config.scenes[model_num]
-        except:
-            return None
 
     def get_usernamespace(self):
         """Namespace only used with #python blocks which are deprecated."""
         tmp = {
             "number_model_runs": sim_config.model_end,
-            "current_model_run": model_num + 1,
+            "current_model_run": self.model_num + 1,
             "inputfile": sim_config.input_file_path.resolve(),
         }
         return dict(**sim_config.em_consts, **tmp)
@@ -194,6 +189,14 @@ class SimulationConfig:
     N.B. A simulation can consist of multiple models.
     """
 
+    # TODO: Make this an enum
+    em_consts = {
+        "c": c,  # Speed of light in free space (m/s)
+        "e0": e0,  # Permittivity of free space (F/m)
+        "m0": m0,  # Permeability of free space (H/m)
+        "z0": np.sqrt(m0 / e0),  # Impedance of free space (Ohms)
+    }
+
     def __init__(self, args):
         """
         Args:
@@ -202,14 +205,38 @@ class SimulationConfig:
 
         self.args = args
 
-        if self.args.taskfarm and self.args.geometry_fixed:
+        self.geometry_fixed: bool = args.geometry_fixed
+        self.geometry_only: bool = args.geometry_only
+        self.gpu: Union[List[str], bool] = args.gpu
+        self.mpi: List[int] = args.mpi
+        self.number_of_models: int = args.n
+        self.opencl: Union[List[str], bool] = args.opencl
+        self.output_file_path: str = args.outputfile
+        self.taskfarm: bool = args.taskfarm
+        self.write_processed_input_file: bool = (
+            args.write_processed
+        )  # For depreciated Python blocks
+
+        if self.taskfarm and self.geometry_fixed:
             logger.exception("The geometry fixed option cannot be used with MPI taskfarm.")
             raise ValueError
 
-        if self.args.gpu and self.args.opencl:
+        if self.gpu and self.opencl:
             logger.exception("You cannot use both CUDA and OpenCl simultaneously.")
             raise ValueError
 
+        if self.mpi and self.args.subgrid:
+            logger.exception("You cannot use subgrids with MPI.")
+            raise ValueError
+
+        # Each model in a simulation is given a unique number when the instance of ModelConfig is created
+        self.current_model = 0
+
+        # Instances of ModelConfig that hold model configuration parameters.
+        # TODO: Consider if this would be better as a dictionary.
+        # Or maybe a non fixed length list (i.e. append each config)
+        self.model_configs: List[Optional[ModelConfig]] = [None] * self.number_of_models
+
         # General settings for the simulation
         #   solver: cpu, cuda, opencl.
         #   precision: data type for electromagnetic field output (single/double).
@@ -217,25 +244,25 @@ class SimulationConfig:
         #                   progressbars when logging level is greater than
         #                   info (20)
 
-        self.general = {"solver": "cpu", "precision": "single", "progressbars": args.log_level <= 20}
+        self.general = {
-
+            "solver": "cpu",
-        self.em_consts = {
+            "precision": "single",
-            "c": c,  # Speed of light in free space (m/s)
+            "progressbars": args.log_level <= 20,
-            "e0": e0,  # Permittivity of free space (F/m)
-            "m0": m0,  # Permeability of free space (H/m)
-            "z0": np.sqrt(m0 / e0),  # Impedance of free space (Ohms)
         }
 
         # Store information about host machine
         self.hostinfo = get_host_info()
 
         # CUDA
-        if self.args.gpu is not None:
+        if self.gpu is not None:
             self.general["solver"] = "cuda"
             # Both single and double precision are possible on GPUs, but single
             # provides best performance.
             self.general["precision"] = "single"
-            self.devices = {"devs": [], "nvcc_opts": None}  # pycuda device objects; nvcc compiler options
+            self.devices = {
+                "devs": [],
+                "nvcc_opts": None,
+            }  # pycuda device objects; nvcc compiler options
             # Suppress nvcc warnings on Microsoft Windows
             if sys.platform == "win32":
                 self.devices["nvcc_opts"] = ["-w"]
@@ -244,14 +271,17 @@ class SimulationConfig:
             self.devices["devs"] = detect_cuda_gpus()
 
         # OpenCL
-        if self.args.opencl is not None:
+        if self.opencl is not None:
             self.general["solver"] = "opencl"
             self.general["precision"] = "single"
-            self.devices = {"devs": [], "compiler_opts": None}  # pyopencl device device(s); compiler options
+            self.devices = {
+                "devs": [],
+                "compiler_opts": None,
+            }  # pyopencl device device(s); compiler options
 
             # Suppress CompilerWarning (sub-class of UserWarning)
             warnings.filterwarnings("ignore", category=UserWarning)
-            
+
             # Suppress unused variable warnings on gcc
             # if sys.platform != 'win32': self.devices['compiler_opts'] = ['-w']
 
@@ -267,43 +297,28 @@ class SimulationConfig:
                 self.general["subgrid"] and self.general["solver"] == "opencl"
             ):
                 logger.exception(
-                    "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids."
+                    "You cannot currently use CUDA or OpenCL-based solvers with models that contain sub-grids."
                 )
                 raise ValueError
         else:
             self.general["subgrid"] = False
 
+        self.autotranslate_subgrid_coordinates = True
+        if hasattr(args, "autotranslate"):
+            self.autotranslate_subgrid_coordinates: bool = args.autotranslate
+
         # Scenes parameter may not exist if user enters via CLI
-        try:
+        self.scenes: List[Optional[Scene]]
-            self.scenes = args.scenes if args.scenes is not None else []
+        if hasattr(args, "scenes") and args.scenes is not None:
-        except AttributeError:
+            self.scenes = args.scenes
-            self.scenes = []
+        else:
+            self.scenes = [None] * self.number_of_models
 
         # Set more complex parameters
         self._set_precision()
         self._set_input_file_path()
         self._set_model_start_end()
 
-    def set_model_device(self, deviceID):
-        """Specify pycuda/pyopencl object for model.
-
-        Args:
-            deviceID: int of requested deviceID of compute device.
-
-        Returns:
-            dev: requested pycuda/pyopencl device object.
-        """
-
-        found = False
-        for ID, dev in self.devices["devs"].items():
-            if ID == deviceID:
-                found = True
-                return dev
-
-        if not found:
-            logger.exception(f"Compute device with device ID {deviceID} does " "not exist.")
-            raise ValueError
-
     def _set_precision(self):
         """Data type (precision) for electromagnetic field output.
 
@@ -342,6 +357,15 @@ class SimulationConfig:
             elif self.general["solver"] == "opencl":
                 self.dtypes["C_complex"] = "cdouble"
 
+    def _set_input_file_path(self):
+        """Sets input file path for CLI or API."""
+        # API
+        if self.args.inputfile is None:
+            self.input_file_path = Path(self.args.outputfile)
+        # API/CLI
+        else:
+            self.input_file_path = Path(self.args.inputfile)
+
     def _set_model_start_end(self):
         """Sets range for number of models to run (internally 0 index)."""
         if self.args.i:
@@ -354,11 +378,94 @@ class SimulationConfig:
         self.model_start = modelstart
         self.model_end = modelend
 
-    def _set_input_file_path(self):
+    def get_model_device(self, deviceID):
-        """Sets input file path for CLI or API."""
+        """Specify pycuda/pyopencl object for model.
-        # API
+
-        if self.args.inputfile is None:
+        Args:
-            self.input_file_path = Path(self.args.outputfile)
+            deviceID: int of requested deviceID of compute device.
-        # API/CLI
+
-        else:
+        Returns:
-            self.input_file_path = Path(self.args.inputfile)
+            dev: requested pycuda/pyopencl device object.
+        """
+
+        found = False
+        for ID, dev in self.devices["devs"].items():
+            if ID == deviceID:
+                found = True
+                return dev
+
+        if not found:
+            logger.exception(f"Compute device with device ID {deviceID} does " "not exist.")
+            raise ValueError
+
+    def get_model_config(self, model_num: Optional[int] = None) -> ModelConfig:
+        """Return ModelConfig instance for specific model.
+
+        Args:
+            model_num: number of the model. If None, returns the config for the current model
+
+        Returns:
+            model_config: requested model config
+        """
+        if model_num is None:
+            model_num = self.current_model
+
+        model_config = self.model_configs[model_num]
+        if model_config is None:
+            logger.exception(f"Cannot get ModelConfig for model {model_num}. It has not been set.")
+            raise ValueError
+
+        return model_config
+
+    def set_model_config(self, model_config: ModelConfig, model_num: Optional[int] = None) -> None:
+        """Set ModelConfig instace for specific model.
+
+        Args:
+            model_num: number of the model. If None, sets the config for the current model
+        """
+        if model_num is None:
+            model_num = self.current_model
+
+        self.model_configs[model_num] = model_config
+
+    def set_current_model(self, model_num: int) -> None:
+        """Set the current model by it's unique identifier
+
+        Args:
+            model_num: unique identifier for the current model
+        """
+        self.current_model = model_num
+
+    def get_scene(self, model_num: Optional[int] = None) -> Optional[Scene]:
+        """Return Scene instance for specific model.
+
+        Args:
+            model_num: number of the model. If None, returns the scene for the current model
+
+        Returns:
+            scene: requested scene
+        """
+        if model_num is None:
+            model_num = self.current_model
+
+        return self.scenes[model_num]
+
+    def set_scene(self, scene: Scene, model_num: Optional[int] = None) -> None:
+        """Set Scene instace for specific model.
+
+        Args:
+            model_num: number of the model. If None, sets the scene for the current model
+        """
+        if model_num is None:
+            model_num = self.current_model
+
+        self.scenes[model_num] = scene
+
+
+# Single instance of SimConfig to hold simulation configuration parameters.
+sim_config: SimulationConfig = None
+
+
+def get_model_config() -> ModelConfig:
+    """Return ModelConfig instance for specific model."""
+    return sim_config.get_model_config()

gprMax/contexts.py

@@ -20,17 +20,22 @@ import datetime
 import gc
 import logging
 import sys
+from typing import Any, Dict, List, Optional
 
 import humanize
 from colorama import Fore, Style, init
 
+from gprMax.hash_cmds_file import parse_hash_commands
+from gprMax.scene import Scene
+
 init()
 
 import gprMax.config as config
+from gprMax.config import ModelConfig
 
 from ._version import __version__, codename
-from .model_build_run import ModelBuildRun
+from .model import Model
-from .solvers import create_G, create_solver
+from .solvers import create_solver
 from .utilities.host_info import print_cuda_info, print_host_info, print_opencl_info
 from .utilities.utilities import get_terminal_width, logo, timer
 
@@ -38,24 +43,24 @@ logger = logging.getLogger(__name__)
 
 
 class Context:
-    """Standard context - models are run one after another and each model
+    """Standard context for building and running models.
-    can exploit parallelisation using either OpenMP (CPU), CUDA (GPU), or
+
-    OpenCL (CPU/GPU).
+    Models are run one after another and each model can exploit
+    parallelisation using either OpenMP (CPU), CUDA (GPU), or OpenCL
+    (CPU/GPU).
     """
 
     def __init__(self):
         self.model_range = range(config.sim_config.model_start, config.sim_config.model_end)
-        self.tsimend = None
+        self.sim_start_time = 0
-        self.tsimstart = None
+        self.sim_end_time = 0
 
-    def run(self):
+    def _start_simulation(self) -> None:
-        """Run the simulation in the correct context.
+        """Run pre-simulation steps
 
-        Returns:
+        Start simulation timer. Output copyright notice and host info.
-            results: dict that can contain useful results/data from simulation.
         """
-
+        self.sim_start_time = timer()
-        self.tsimstart = timer()
         self.print_logo_copyright()
         print_host_info(config.sim_config.hostinfo)
         if config.sim_config.general["solver"] == "cuda":
@@ -63,56 +68,128 @@ class Context:
         elif config.sim_config.general["solver"] == "opencl":
             print_opencl_info(config.sim_config.devices["devs"])
 
-        # Clear list of model configs. It can be retained when gprMax is
+    def _end_simulation(self) -> None:
-        # called in a loop, and want to avoid this.
+        """Run post-simulation steps
-        config.model_configs = []
+
+        Stop simulation timer. Output timing information.
+        """
+        self.sim_end_time = timer()
+        self.print_sim_time_taken()
+
+    def run(self) -> Dict:
+        """Run the simulation in the correct context.
+
+        Returns:
+            results: dict that can contain useful results/data from simulation.
+        """
+
+        self._start_simulation()
 
         for i in self.model_range:
-            config.model_num = i
+            self._run_model(i)
-            model_config = config.ModelConfig()
-            config.model_configs.append(model_config)
 
-            # Always create a grid for the first model. The next model to run
+        self._end_simulation()
-            # only gets a new grid if the geometry is not re-used.
-            if i != 0 and config.sim_config.args.geometry_fixed:
-                config.get_model_config().reuse_geometry = True
-            else:
-                G = create_G()
-
-            model = ModelBuildRun(G)
-            model.build()
-
-            if not config.sim_config.args.geometry_only:
-                solver = create_solver(G)
-                model.solve(solver)
-                del solver, model
-
-            if not config.sim_config.args.geometry_fixed:
-                # Manual garbage collection required to stop memory leak on GPUs
-                # when using pycuda
-                del G
-
-            gc.collect()
-
-        self.tsimend = timer()
-        self.print_sim_time_taken()
 
         return {}
 
-    def print_logo_copyright(self):
+    def _run_model(self, model_num: int) -> None:
+        """Process for running a single model.
+
+        Args:
+            model_num: index of model to be run
+        """
+
+        config.sim_config.set_current_model(model_num)
+        model_config = self._create_model_config(model_num)
+        config.sim_config.set_model_config(model_config)
+
+        if not model_config.reuse_geometry():
+            scene = self._get_scene(model_num)
+            model = self._create_model()
+            scene.create_internal_objects(model.G)
+
+        model.build()
+
+        if not config.sim_config.geometry_only:
+            solver = create_solver(model.G)
+            model.solve(solver)
+            del solver
+
+        if not config.sim_config.geometry_fixed:
+            # Manual garbage collection required to stop memory leak on GPUs
+            # when using pycuda
+            del model.G, model
+
+        gc.collect()
+
+    def _create_model_config(self, model_num: int) -> ModelConfig:
+        """Create model config and save to global config."""
+        return ModelConfig(model_num)
+
+    def _get_scene(self, model_num: int) -> Scene:
+        # API for multiple scenes / model runs
+        scene = config.sim_config.get_scene(model_num)
+
+        # If there is no scene, process the hash commands
+        if scene is None:
+            scene = Scene()
+            config.sim_config.set_scene(scene, model_num)
+            # Parse the input file into user objects and add them to the scene
+            scene = parse_hash_commands(scene)
+
+        return scene
+
+    def _create_model(self) -> Model:
+        return Model()
+
+    def print_logo_copyright(self) -> None:
         """Prints gprMax logo, version, and copyright/licencing information."""
         logo_copyright = logo(f"{__version__} ({codename})")
         logger.basic(logo_copyright)
 
-    def print_sim_time_taken(self):
+    def print_sim_time_taken(self) -> None:
         """Prints the total simulation time based on context."""
         s = (
             f"\n=== Simulation completed in "
-            f"{humanize.precisedelta(datetime.timedelta(seconds=self.tsimend - self.tsimstart), format='%0.4f')}"
+            f"{humanize.precisedelta(datetime.timedelta(seconds=self.sim_end_time - self.sim_start_time), format='%0.4f')}"
         )
         logger.basic(f"{s} {'=' * (get_terminal_width() - 1 - len(s))}\n")
 
 
+class MPIContext(Context):
+    def __init__(self, x_dim: int, y_dim: int, z_dim: int):
+        super().__init__()
+        from mpi4py import MPI
+
+        self.comm = MPI.COMM_WORLD
+        self.rank = self.comm.rank
+
+        if self.rank >= x_dim * y_dim * z_dim:
+            logger.warn(
+                (
+                    f"Rank {self.rank}: Only {x_dim * y_dim * z_dim} MPI ranks required for the"
+                    " dimensions specified. Either increase your MPI dimension size, or request"
+                    " fewer MPI tasks."
+                )
+            )
+            self.x = -1
+            self.y = -1
+            self.z = -1
+        else:
+            self.x = self.rank % x_dim
+            self.y = (self.rank // x_dim) % y_dim
+            self.z = (self.rank // (x_dim * y_dim)) % z_dim
+
+    def run(self) -> Dict:
+        print(f"I am rank {self.rank} and I will run at grid position {self.x}, {self.y}, {self.z}")
+
+        if self.rank == 0:
+            print("Rank 0 is running the simulation")
+            return super().run()
+        else:
+            return {}
+
+
 class TaskfarmContext(Context):
     """Mixed mode MPI/OpenMP/CUDA context - MPI task farm is used to distribute
     models, and each model parallelised using either OpenMP (CPU),
@@ -129,17 +206,12 @@ class TaskfarmContext(Context):
         self.rank = self.comm.rank
         self.TaskfarmExecutor = TaskfarmExecutor
 
-    def _run_model(self, **work):
+    def _create_model_config(self, model_num: int) -> ModelConfig:
-        """Process for running a single model.
+        """Create model config and save to global config.
 
-        Args:
+        Set device in model config according to MPI rank.
-            work: dict of any additional information that is passed to MPI
-                    workers. By default only model number (i) is used.
         """
-
+        model_config = super()._create_model_config(model_num)
-        # Create configuration for model
-        config.model_num = work["i"]
-        model_config = config.ModelConfig()
         # Set GPU deviceID according to worker rank
         if config.sim_config.general["solver"] == "cuda":
             model_config.device = {
@@ -147,23 +219,18 @@ class TaskfarmContext(Context):
                 "deviceID": self.rank - 1,
                 "snapsgpu2cpu": False,
             }
-        config.model_configs = model_config
+        return model_config
 
-        G = create_G()
+    def _run_model(self, **work) -> None:
-        model = ModelBuildRun(G)
+        """Process for running a single model.
-        model.build()
 
-        if not config.sim_config.args.geometry_only:
+        Args:
-            solver = create_solver(G)
+            work: dict of any additional information that is passed to MPI
-            model.solve(solver)
+                workers. By default only model number (i) is used.
-            del solver, model
+        """
+        return super()._run_model(work["i"])
 
-        # Manual garbage collection required to stop memory leak on GPUs when
+    def run(self) -> Optional[List[Optional[Dict]]]:
-        # using pycuda
-        del G
-        gc.collect()
-
-    def run(self):
         """Specialise how the models are run.
 
         Returns:
@@ -171,16 +238,12 @@ class TaskfarmContext(Context):
         """
 
         if self.rank == 0:
-            self.tsimstart = timer()
+            self._start_simulation()
-            self.print_logo_copyright()
-            print_host_info(config.sim_config.hostinfo)
-            if config.sim_config.general["solver"] == "cuda":
-                print_cuda_info(config.sim_config.devices["devs"])
-            elif config.sim_config.general["solver"] == "opencl":
-                print_opencl_info(config.sim_config.devices["devs"])
 
             s = f"\n--- Input file: {config.sim_config.input_file_path}"
-            logger.basic(Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL)
+            logger.basic(
+                Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
+            )
 
             sys.stdout.flush()
 
@@ -210,6 +273,5 @@ class TaskfarmContext(Context):
         executor.join()
 
         if executor.is_master():
-            self.tsimend = timer()
+            self._end_simulation()
-            self.print_sim_time_taken()
             return results

gprMax/gprMax.py

@@ -15,12 +15,11 @@
 #
 # 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 gprMax.config as config
 
-from .contexts import Context, TaskfarmContext
+from .contexts import Context, MPIContext, TaskfarmContext
 from .utilities.logging import logging_config
 
 # Arguments (used for API) and their default values (used for API and CLI)
@@ -31,6 +30,7 @@ args_defaults = {
     "n": 1,
     "i": None,
     "taskfarm": False,
+    "mpi": None,
     "gpu": None,
     "opencl": None,
     "subgrid": False,
@@ -44,34 +44,53 @@ args_defaults = {
 
 # Argument help messages (used for CLI argparse)
 help_msg = {
-    "scenes": "(list, req): Scenes to run the model. Multiple scene objects "
+    "scenes": (
-    "can given in order to run multiple simulation runs. Each scene "
+        "(list, req): Scenes to run the model. Multiple scene objects can given in order to run"
-    "must contain the essential simulation objects",
+        " multiple simulation runs. Each scene must contain the essential simulation objects"
-    "inputfile": "(str, opt): Input file path. Can also run simulation by " "providing an input file.",
+    ),
+    "inputfile": "(str, opt): Input file path. Can also run simulation by providing an input file.",
     "outputfile": "(str, req): File path to the output data file.",
     "n": "(int, req): Number of required simulation runs.",
-    "i": "(int, opt): Model number to start/restart simulation from. It would "
+    "i": (
-    "typically be used to restart a series of models from a specific "
+        "(int, opt): Model number to start/restart simulation from. It would typically be used to"
-    "model number, with the n argument, e.g. to restart from A-scan 45 "
+        " restart a series of models from a specific model number, with the n argument, e.g. to"
-    "when creating a B-scan with 60 traces.",
+        " restart from A-scan 45 when creating a B-scan with 60 traces."
-    "taskfarm": "(bool, opt): Flag to use Message Passing Interface (MPI) task farm. "
+    ),
-    "This option is most usefully combined with n to allow individual "
+    "taskfarm": (
-    "models to be farmed out using a MPI task farm, e.g. to create a "
+        "(bool, opt): Flag to use Message Passing Interface (MPI) task farm. This option is most"
-    "B-scan with 60 traces and use MPI to farm out each trace. For "
+        " usefully combined with n to allow individual models to be farmed out using a MPI task"
-    "further details see the performance section of the User Guide.",
+        " farm, e.g. to create a B-scan with 60 traces and use MPI to farm out each trace. For"
-    "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA GPU " "device ID(s) for specific GPU card(s).",
+        " further details see the performance section of the User Guide."
-    "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL device " "ID(s) for specific compute device(s).",
+    ),
+    "mpi": (
+        "(list, opt): Flag to use Message Passing Interface (MPI) to divide the model between MPI"
+        " ranks. Three integers should be provided to define the number of MPI processes (min 1) in"
+        " the x, y, and z dimensions."
+    ),
+    "gpu": (
+        "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA GPU device ID(s) for specific"
+        " GPU card(s)."
+    ),
+    "opencl": (
+        "(list/bool, opt): Flag to use OpenCL or list of OpenCL device ID(s) for specific compute"
+        " device(s)."
+    ),
     "subgrid": "(bool, opt): Flag to use sub-gridding.",
-    "autotranslate": "(bool, opt): For sub-gridding - auto translate objects "
+    "autotranslate": (
-    "with main grid coordinates to their equivalent local "
+        "(bool, opt): For sub-gridding - auto translate objects with main grid coordinates to their"
-    "grid coordinate within the subgrid. If this option is "
+        " equivalent local grid coordinate within the subgrid. If this option is off users must"
-    "off users must specify sub-grid object point within the "
+        " specify sub-grid object point within the global subgrid space."
-    "global subgrid space.",
+    ),
-    "geometry_only": "(bool, opt): Build a model and produce any geometry " "views but do not run the simulation.",
+    "geometry_only": (
-    "geometry_fixed": "(bool, opt): Run a series of models where the geometry " "does not change between models.",
+        "(bool, opt): Build a model and produce any geometry views but do not run the simulation."
-    "write_processed": "(bool, opt): Writes another input file after any "
+    ),
-    "Python code (#python blocks) and in the original input "
+    "geometry_fixed": (
-    "file has been processed.",
+        "(bool, opt): Run a series of models where the geometry does not change between models."
+    ),
+    "write_processed": (
+        "(bool, opt): Writes another input file after any Python code (#python blocks) and in the"
+        " original input file has been processed."
+    ),
     "log_level": "(int, opt): Level of logging to use.",
     "log_file": "(bool, opt): Write logging information to file.",
 }
@@ -84,6 +103,7 @@ def run(
     n=args_defaults["n"],
     i=args_defaults["i"],
     taskfarm=args_defaults["taskfarm"],
+    mpi=args_defaults["mpi"],
     gpu=args_defaults["gpu"],
     opencl=args_defaults["opencl"],
     subgrid=args_defaults["subgrid"],
@@ -94,44 +114,51 @@ def run(
     log_level=args_defaults["log_level"],
     log_file=args_defaults["log_file"],
 ):
-    """Entry point for application programming interface (API). Runs the
+    """Entry point for application programming interface (API).
-        simulation for the given list of scenes.
+
+    Runs the simulation for the given list of scenes.
 
     Args:
-        scenes: list of the scenes to run the model. Multiple scene objects can
+        scenes: list of the scenes to run the model. Multiple scene
-                be given in order to run multiple simulation runs. Each scene
+            objects can be given in order to run multiple simulation
-                must contain the essential simulation objects.
+            runs. Each scene must contain the essential simulation
-        inputfile: optional string for input file path. Can also run simulation
+            objects.
-                    by providing an input file.
+        inputfile: optional string for input file path. Can also run
+            simulation by providing an input file.
         outputfile: string for file path to the output data file
         n: optional int for number of required simulation runs.
-        i: optional int for model number to start/restart simulation from.
+        i: optional int for model number to start/restart simulation
-            It would typically be used to restart a series of models from a
+            from. It would typically be used to restart a series of
-            specific model number, with the n argument, e.g. to restart from
+            models from a specific model number, with the n argument,
-            A-scan 45 when creating a B-scan with 60 traces.
+            e.g. to restart from A-scan 45 when creating a B-scan with
-        taskfarm: optional boolean flag to use Message Passing Interface (MPI) task
+            60 traces.
-                farm. This option is most usefully combined with n to allow
+        taskfarm: optional boolean flag to use Message Passing Interface
-                individual models to be farmed out using a MPI task farm,
+            (MPI) task farm. This option is most usefully combined with
-                e.g. to create a B-scan with 60 traces and use MPI to farm out
+            n to allow individual models to be farmed out using a MPI
-                each trace. For further details see the performance section of
+            task farm, e.g. to create a B-scan with 60 traces and use
-                the User Guide
+            MPI to farm out each trace. For further details see the
-        gpu: optional list/boolean to use NVIDIA GPU or list of NVIDIA GPU device
+            performance section of the User Guide.
-                ID(s) for specific GPU card(s).
+        mpi: optional flag to use Message Passing Interface (MPI) to
-        opencl: optional list/boolean to use OpenCL or list of OpenCL device ID(s)
+            divide the model between MPI ranks. Three integers should be
-                for specific compute device(s).
+            provided to define the number of MPI processes (min 1) in
+            the x, y, and z dimensions.
+        gpu: optional list/boolean to use NVIDIA GPU or list of NVIDIA
+            GPU device ID(s) for specific GPU card(s).
+        opencl: optional list/boolean to use OpenCL or list of OpenCL
+            device ID(s) for specific compute device(s).
         subgrid: optional boolean to use sub-gridding.
-        autotranslate: optional boolean for sub-gridding to auto translate
+        autotranslate: optional boolean for sub-gridding to auto
-                        objects with main grid coordinates to their equivalent
+            translate objects with main grid coordinates to their
-                        local grid coordinate within the subgrid. If this option
+            equivalent local grid coordinate within the subgrid. If this
-                        is off users must specify sub-grid object point within
+            option is off users must specify sub-grid object point
-                        the global subgrid space.
+            within the global subgrid space.
         geometry_only: optional boolean to build a model and produce any
-                        geometry views but do not run the simulation.
+            geometry views but do not run the simulation.
-        geometry_fixed: optional boolean to run a series of models where the
+        geometry_fixed: optional boolean to run a series of models where
-                        geometry does not change between models.
+            the geometry does not change between models.
-        write_processed: optional boolean to write another input file after any
+        write_processed: optional boolean to write another input file
-                            #python blocks (which are deprecated) in the
+            after any #python blocks (which are deprecated) in the
-                            original input file has been processed.
+            original input file has been processed.
         log_level: optional int for level of logging to use.
         log_file: optional boolean to write logging information to file.
     """
@@ -144,6 +171,7 @@ def run(
             "n": n,
             "i": i,
             "taskfarm": taskfarm,
+            "mpi": mpi,
             "gpu": gpu,
             "opencl": opencl,
             "subgrid": subgrid,
@@ -163,16 +191,34 @@ def cli():
     """Entry point for command line interface (CLI)."""
 
     # Parse command line arguments
-    parser = argparse.ArgumentParser(prog="gprMax", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    parser = argparse.ArgumentParser(
+        prog="gprMax", formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
     parser.add_argument("inputfile", help=help_msg["inputfile"])
     parser.add_argument("-outputfile", "-o", help=help_msg["outputfile"])
     parser.add_argument("-n", default=args_defaults["n"], type=int, help=help_msg["n"])
     parser.add_argument("-i", type=int, help=help_msg["i"])
-    parser.add_argument("-taskfarm", action="store_true", default=args_defaults["taskfarm"], help=help_msg["taskfarm"])
+    parser.add_argument(
+        "-taskfarm",
+        action="store_true",
+        default=args_defaults["taskfarm"],
+        help=help_msg["taskfarm"],
+    )
+    parser.add_argument(
+        "-mpi",
+        type=int,
+        action="store",
+        nargs=3,
+        default=args_defaults["mpi"],
+        help=help_msg["mpi"],
+    )
     parser.add_argument("-gpu", type=int, action="append", nargs="*", help=help_msg["gpu"])
     parser.add_argument("-opencl", type=int, action="append", nargs="*", help=help_msg["opencl"])
     parser.add_argument(
-        "--geometry-only", action="store_true", default=args_defaults["geometry_only"], help=help_msg["geometry_only"]
+        "--geometry-only",
+        action="store_true",
+        default=args_defaults["geometry_only"],
+        help=help_msg["geometry_only"],
     )
     parser.add_argument(
         "--geometry-fixed",
@@ -186,8 +232,15 @@ def cli():
         default=args_defaults["write_processed"],
         help=help_msg["write_processed"],
     )
-    parser.add_argument("--log-level", type=int, default=args_defaults["log_level"], help=help_msg["log_level"])
+    parser.add_argument(
-    parser.add_argument("--log-file", action="store_true", default=args_defaults["log_file"], help=help_msg["log_file"])
+        "--log-level", type=int, default=args_defaults["log_level"], help=help_msg["log_level"]
+    )
+    parser.add_argument(
+        "--log-file",
+        action="store_true",
+        default=args_defaults["log_file"],
+        help=help_msg["log_file"],
+    )
     args = parser.parse_args()
 
     results = run_main(args)
@@ -202,8 +255,9 @@ def run_main(args):
         args: namespace with arguments from either API or CLI.
 
     Returns:
-        results: dict that can contain useful results/data from simulation.
+        results: dict that can contain useful results/data from
-                    Enables these to be propagated to calling script.
+            simulation. Enables these to be propagated to calling
+            script.
     """
 
     results = {}
@@ -213,6 +267,12 @@ def run_main(args):
     # MPI taskfarm running with (OpenMP/CUDA/OpenCL)
     if config.sim_config.args.taskfarm:
         context = TaskfarmContext()
+    # MPI running to divide model between ranks
+    elif config.sim_config.args.mpi is not None:
+        x = config.sim_config.args.mpi[0]
+        y = config.sim_config.args.mpi[1]
+        z = config.sim_config.args.mpi[2]
+        context = MPIContext(x, y, z)
     # Standard running (OpenMP/CUDA/OpenCL)
     else:
         context = Context()

gprMax/grid/cuda_grid.py

@@ -0,0 +1,67 @@
+# Copyright (C) 2015-2024: 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/>.
+
+from importlib import import_module
+
+import numpy as np
+
+from gprMax.grid.fdtd_grid import FDTDGrid
+
+
+class CUDAGrid(FDTDGrid):
+    """Additional grid methods for solving on GPU using CUDA."""
+
+    def __init__(self):
+        super().__init__()
+
+        self.gpuarray = import_module("pycuda.gpuarray")
+
+        # Threads per block - used for main electric/magnetic field updates
+        self.tpb = (128, 1, 1)
+        # Blocks per grid - used for main electric/magnetic field updates
+        self.bpg = None
+
+    def set_blocks_per_grid(self):
+        """Set the blocks per grid size used for updating the electric and
+        magnetic field arrays on a GPU.
+        """
+
+        self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), 1, 1)
+
+    def htod_geometry_arrays(self):
+        """Initialise an array for cell edge IDs (ID) on compute device."""
+
+        self.ID_dev = self.gpuarray.to_gpu(self.ID)
+
+    def htod_field_arrays(self):
+        """Initialise field arrays on compute device."""
+
+        self.Ex_dev = self.gpuarray.to_gpu(self.Ex)
+        self.Ey_dev = self.gpuarray.to_gpu(self.Ey)
+        self.Ez_dev = self.gpuarray.to_gpu(self.Ez)
+        self.Hx_dev = self.gpuarray.to_gpu(self.Hx)
+        self.Hy_dev = self.gpuarray.to_gpu(self.Hy)
+        self.Hz_dev = self.gpuarray.to_gpu(self.Hz)
+
+    def htod_dispersive_arrays(self):
+        """Initialise dispersive material coefficient arrays on compute device."""
+
+        self.updatecoeffsdispersive_dev = self.gpuarray.to_gpu(self.updatecoeffsdispersive)
+        self.Tx_dev = self.gpuarray.to_gpu(self.Tx)
+        self.Ty_dev = self.gpuarray.to_gpu(self.Ty)
+        self.Tz_dev = self.gpuarray.to_gpu(self.Tz)

gprMax/grid/fdtd_grid.py

@@ -16,18 +16,32 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-import decimal as d
+import decimal
+import itertools
+import logging
+import sys
 from collections import OrderedDict
-from importlib import import_module
+from typing import Iterable, List, Union
 
+import humanize
 import numpy as np
+from terminaltables import SingleTable
+from tqdm import tqdm
 
-import gprMax.config as config
+from gprMax import config
+from gprMax.cython.yee_cell_build import build_electric_components, build_magnetic_components
 
-from .pml import PML
+# from gprMax.geometry_outputs import GeometryObjects, GeometryView
-from .utilities.utilities import fft_power, round_value
+from gprMax.materials import Material, process_materials
+from gprMax.pml import CFS, PML, build_pml, print_pml_info
+from gprMax.receivers import Rx
+from gprMax.sources import HertzianDipole, MagneticDipole, Source, VoltageSource
 
-np.seterr(invalid="raise")
+# from gprMax.subgrids.grid import SubGridBaseGrid
+from gprMax.utilities.host_info import mem_check_build_all, mem_check_run_all
+from gprMax.utilities.utilities import fft_power, get_terminal_width, round_value
+
+logger = logging.getLogger(__name__)
 
 
 class FDTDGrid:
@@ -64,23 +78,218 @@ class FDTDGrid:
         # corrections will be different.
         self.pmls["thickness"] = OrderedDict((key, 10) for key in PML.boundaryIDs)
 
-        self.materials = []
+        # TODO: Add type information.
+        # Currently importing GeometryObjects, GeometryView, and
+        # SubGridBaseGrid cause cyclic dependencies
+        self.materials: List[Material] = []
         self.mixingmodels = []
         self.averagevolumeobjects = True
         self.fractalvolumes = []
         self.geometryviews = []
         self.geometryobjectswrite = []
         self.waveforms = []
-        self.voltagesources = []
+        self.voltagesources: List[VoltageSource] = []
-        self.hertziandipoles = []
+        self.hertziandipoles: List[HertzianDipole] = []
-        self.magneticdipoles = []
+        self.magneticdipoles: List[MagneticDipole] = []
         self.transmissionlines = []
-        self.rxs = []
+        self.rxs: List[Rx] = []
-        self.srcsteps = [0, 0, 0]
+        self.srcsteps: List[float] = [0, 0, 0]
-        self.rxsteps = [0, 0, 0]
+        self.rxsteps: List[float] = [0, 0, 0]
         self.snapshots = []
         self.subgrids = []
 
+    def build(self) -> None:
+        # Print info on any subgrids
+        for subgrid in self.subgrids:
+            subgrid.print_info()
+
+        # Combine available grids
+        grids = [self] + self.subgrids
+
+        # Check for dispersive materials (and specific type)
+        if config.get_model_config().materials["maxpoles"] != 0:
+            # TODO: This sets materials["drudelorentz"] based only the
+            # last grid/subgrid. Is that correct?
+            for grid in grids:
+                config.get_model_config().materials["drudelorentz"] = any(
+                    [m for m in grid.materials if "drude" in m.type or "lorentz" in m.type]
+                )
+
+            # Set data type if any dispersive materials (must be done before memory checks)
+            config.get_model_config().set_dispersive_material_types()
+
+        # Check memory requirements to build model/scene (different to memory
+        # requirements to run model when FractalVolumes/FractalSurfaces are
+        # used as these can require significant additional memory)
+        total_mem_build, mem_strs_build = mem_check_build_all(grids)
+
+        # Check memory requirements to run model
+        total_mem_run, mem_strs_run = mem_check_run_all(grids)
+
+        if total_mem_build > total_mem_run:
+            logger.info(
+                f'\nMemory required (estimated): {" + ".join(mem_strs_build)} + '
+                f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
+                f"overhead = {humanize.naturalsize(total_mem_build)}"
+            )
+        else:
+            logger.info(
+                f'\nMemory required (estimated): {" + ".join(mem_strs_run)} + '
+                f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
+                f"overhead = {humanize.naturalsize(total_mem_run)}"
+            )
+
+        # Build grids
+        for grid in grids:
+            # Set default CFS parameter for PMLs if not user provided
+            if not grid.pmls["cfs"]:
+                grid.pmls["cfs"] = [CFS()]
+            logger.info(print_pml_info(grid))
+            if not all(value == 0 for value in grid.pmls["thickness"].values()):
+                grid._build_pmls()
+            if grid.averagevolumeobjects:
+                grid._build_components()
+            grid._tm_grid_update()
+            grid._update_voltage_source_materials()
+            grid.initialise_field_arrays()
+            grid.initialise_std_update_coeff_arrays()
+            if config.get_model_config().materials["maxpoles"] > 0:
+                grid.initialise_dispersive_arrays()
+                grid.initialise_dispersive_update_coeff_array()
+            grid._build_materials()
+
+            # Check to see if numerical dispersion might be a problem
+            results = dispersion_analysis(grid)
+            if results["error"]:
+                logger.warning(
+                    f"\nNumerical dispersion analysis [{grid.name}] "
+                    f"not carried out as {results['error']}"
+                )
+            elif results["N"] < config.get_model_config().numdispersion["mingridsampling"]:
+                logger.exception(
+                    f"\nNon-physical wave propagation in [{grid.name}] "
+                    f"detected. Material '{results['material'].ID}' "
+                    f"has wavelength sampled by {results['N']} cells, "
+                    f"less than required minimum for physical wave "
+                    f"propagation. Maximum significant frequency "
+                    f"estimated as {results['maxfreq']:g}Hz"
+                )
+                raise ValueError
+            elif (
+                results["deltavp"]
+                and np.abs(results["deltavp"])
+                > config.get_model_config().numdispersion["maxnumericaldisp"]
+            ):
+                logger.warning(
+                    f"\n[{grid.name}] has potentially significant "
+                    f"numerical dispersion. Estimated largest physical "
+                    f"phase-velocity error is {results['deltavp']:.2f}% "
+                    f"in material '{results['material'].ID}' whose "
+                    f"wavelength sampled by {results['N']} cells. "
+                    f"Maximum significant frequency estimated as "
+                    f"{results['maxfreq']:g}Hz"
+                )
+            elif results["deltavp"]:
+                logger.info(
+                    f"\nNumerical dispersion analysis [{grid.name}]: "
+                    f"estimated largest physical phase-velocity error is "
+                    f"{results['deltavp']:.2f}% in material '{results['material'].ID}' "
+                    f"whose wavelength sampled by {results['N']} cells. "
+                    f"Maximum significant frequency estimated as "
+                    f"{results['maxfreq']:g}Hz"
+                )
+
+    def _build_pmls(self) -> None:
+        pbar = tqdm(
+            total=sum(1 for value in self.pmls["thickness"].values() if value > 0),
+            desc=f"Building PML boundaries [{self.name}]",
+            ncols=get_terminal_width() - 1,
+            file=sys.stdout,
+            disable=not config.sim_config.general["progressbars"],
+        )
+        for pml_id, thickness in self.pmls["thickness"].items():
+            if thickness > 0:
+                build_pml(self, pml_id, thickness)
+                pbar.update()
+        pbar.close()
+
+    def _build_components(self) -> None:
+        # Build the model, i.e. set the material properties (ID) for every edge
+        # of every Yee cell
+        logger.info("")
+        pbar = tqdm(
+            total=2,
+            desc=f"Building Yee cells [{self.name}]",
+            ncols=get_terminal_width() - 1,
+            file=sys.stdout,
+            disable=not config.sim_config.general["progressbars"],
+        )
+        build_electric_components(self.solid, self.rigidE, self.ID, self)
+        pbar.update()
+        build_magnetic_components(self.solid, self.rigidH, self.ID, self)
+        pbar.update()
+        pbar.close()
+
+    def _tm_grid_update(self) -> None:
+        if config.get_model_config().mode == "2D TMx":
+            self.tmx()
+        elif config.get_model_config().mode == "2D TMy":
+            self.tmy()
+        elif config.get_model_config().mode == "2D TMz":
+            self.tmz()
+
+    def _update_voltage_source_materials(self):
+        # Process any voltage sources (that have resistance) to create a new
+        # material at the source location
+        for voltagesource in self.voltagesources:
+            voltagesource.create_material(self)
+
+    def _build_materials(self) -> None:
+        # Process complete list of materials - calculate update coefficients,
+        # store in arrays, and build text list of materials/properties
+        materialsdata = process_materials(self)
+        materialstable = SingleTable(materialsdata)
+        materialstable.outer_border = False
+        materialstable.justify_columns[0] = "right"
+
+        logger.info(f"\nMaterials [{self.name}]:")
+        logger.info(materialstable.table)
+
+    def _update_positions(
+        self, items: Iterable[Union[Source, Rx]], step_size: List[float], step_number: int
+    ) -> None:
+        if step_size[0] != 0 or step_size[1] != 0 or step_size[2] != 0:
+            for item in items:
+                if step_number == 0:
+                    if (
+                        item.xcoord + self.srcsteps[0] * config.sim_config.model_end < 0
+                        or item.xcoord + self.srcsteps[0] * config.sim_config.model_end > self.nx
+                        or item.ycoord + self.srcsteps[1] * config.sim_config.model_end < 0
+                        or item.ycoord + self.srcsteps[1] * config.sim_config.model_end > self.ny
+                        or item.zcoord + self.srcsteps[2] * config.sim_config.model_end < 0
+                        or item.zcoord + self.srcsteps[2] * config.sim_config.model_end > self.nz
+                    ):
+                        raise ValueError
+                item.xcoord = item.xcoordorigin + step_number * step_size[0]
+                item.ycoord = item.ycoordorigin + step_number * step_size[1]
+                item.zcoord = item.zcoordorigin + step_number * step_size[2]
+
+    def update_simple_source_positions(self, step: int = 0) -> None:
+        try:
+            self._update_positions(
+                itertools.chain(self.hertziandipoles, self.magneticdipoles), self.srcsteps, step
+            )
+        except ValueError as e:
+            logger.exception("Source(s) will be stepped to a position outside the domain.")
+            raise ValueError from e
+
+    def update_receiver_positions(self, step: int = 0) -> None:
+        try:
+            self._update_positions(self.rxs, self.rxsteps, step)
+        except ValueError as e:
+            logger.exception("Receiver(s) will be stepped to a position outside the domain.")
+            raise ValueError from e
+
     def within_bounds(self, p):
         if p[0] < 0 or p[0] > self.nx:
             raise ValueError("x")
@@ -128,30 +337,67 @@ class FDTDGrid:
 
     def initialise_field_arrays(self):
         """Initialise arrays for the electric and magnetic field components."""
-        self.Ex = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+        self.Ex = np.zeros(
-        self.Ey = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+            (self.nx + 1, self.ny + 1, self.nz + 1),
-        self.Ez = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+            dtype=config.sim_config.dtypes["float_or_double"],
-        self.Hx = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+        )
-        self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+        self.Ey = np.zeros(
-        self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"])
+            (self.nx + 1, self.ny + 1, self.nz + 1),
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
+        self.Ez = np.zeros(
+            (self.nx + 1, self.ny + 1, self.nz + 1),
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
+        self.Hx = np.zeros(
+            (self.nx + 1, self.ny + 1, self.nz + 1),
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
+        self.Hy = np.zeros(
+            (self.nx + 1, self.ny + 1, self.nz + 1),
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
+        self.Hz = np.zeros(
+            (self.nx + 1, self.ny + 1, self.nz + 1),
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
 
     def initialise_std_update_coeff_arrays(self):
         """Initialise arrays for storing update coefficients."""
-        self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=config.sim_config.dtypes["float_or_double"])
+        self.updatecoeffsE = np.zeros(
-        self.updatecoeffsH = np.zeros((len(self.materials), 5), dtype=config.sim_config.dtypes["float_or_double"])
+            (len(self.materials), 5), dtype=config.sim_config.dtypes["float_or_double"]
+        )
+        self.updatecoeffsH = np.zeros(
+            (len(self.materials), 5), dtype=config.sim_config.dtypes["float_or_double"]
+        )
 
     def initialise_dispersive_arrays(self):
         """Initialise field arrays when there are dispersive materials present."""
         self.Tx = np.zeros(
-            (config.get_model_config().materials["maxpoles"], self.nx + 1, self.ny + 1, self.nz + 1),
+            (
+                config.get_model_config().materials["maxpoles"],
+                self.nx + 1,
+                self.ny + 1,
+                self.nz + 1,
+            ),
             dtype=config.get_model_config().materials["dispersivedtype"],
         )
         self.Ty = np.zeros(
-            (config.get_model_config().materials["maxpoles"], self.nx + 1, self.ny + 1, self.nz + 1),
+            (
+                config.get_model_config().materials["maxpoles"],
+                self.nx + 1,
+                self.ny + 1,
+                self.nz + 1,
+            ),
             dtype=config.get_model_config().materials["dispersivedtype"],
         )
         self.Tz = np.zeros(
-            (config.get_model_config().materials["maxpoles"], self.nx + 1, self.ny + 1, self.nz + 1),
+            (
+                config.get_model_config().materials["maxpoles"],
+                self.nx + 1,
+                self.ny + 1,
+                self.nz + 1,
+            ),
             dtype=config.get_model_config().materials["dispersivedtype"],
         )
 
@@ -288,111 +534,29 @@ class FDTDGrid:
     def calculate_dt(self):
         """Calculate time step at the CFL limit."""
         if config.get_model_config().mode == "2D TMx":
-            self.dt = 1 / (config.sim_config.em_consts["c"] * np.sqrt((1 / self.dy**2) + (1 / self.dz**2)))
+            self.dt = 1 / (
+                config.sim_config.em_consts["c"] * np.sqrt((1 / self.dy**2) + (1 / self.dz**2))
+            )
         elif config.get_model_config().mode == "2D TMy":
-            self.dt = 1 / (config.sim_config.em_consts["c"] * np.sqrt((1 / self.dx**2) + (1 / self.dz**2)))
+            self.dt = 1 / (
+                config.sim_config.em_consts["c"] * np.sqrt((1 / self.dx**2) + (1 / self.dz**2))
+            )
         elif config.get_model_config().mode == "2D TMz":
-            self.dt = 1 / (config.sim_config.em_consts["c"] * np.sqrt((1 / self.dx**2) + (1 / self.dy**2)))
+            self.dt = 1 / (
+                config.sim_config.em_consts["c"] * np.sqrt((1 / self.dx**2) + (1 / self.dy**2))
+            )
         else:
             self.dt = 1 / (
-                config.sim_config.em_consts["c"] * np.sqrt((1 / self.dx**2) + (1 / self.dy**2) + (1 / self.dz**2))
+                config.sim_config.em_consts["c"]
+                * np.sqrt((1 / self.dx**2) + (1 / self.dy**2) + (1 / self.dz**2))
             )
 
         # Round down time step to nearest float with precision one less than
         # hardware maximum. Avoids inadvertently exceeding the CFL due to
         # binary representation of floating point number.
-        self.dt = round_value(self.dt, decimalplaces=d.getcontext().prec - 1)
+        self.dt = round_value(self.dt, decimalplaces=decimal.getcontext().prec - 1)
 
 
-class CUDAGrid(FDTDGrid):
-    """Additional grid methods for solving on GPU using CUDA."""
-
-    def __init__(self):
-        super().__init__()
-
-        self.gpuarray = import_module("pycuda.gpuarray")
-
-        # Threads per block - used for main electric/magnetic field updates
-        self.tpb = (128, 1, 1)
-        # Blocks per grid - used for main electric/magnetic field updates
-        self.bpg = None
-
-    def set_blocks_per_grid(self):
-        """Set the blocks per grid size used for updating the electric and
-            magnetic field arrays on a GPU.
-        """
-
-        self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), 1, 1)
-
-    def htod_geometry_arrays(self):
-        """Initialise an array for cell edge IDs (ID) on compute device."""
-
-        self.ID_dev = self.gpuarray.to_gpu(self.ID)
-
-    def htod_field_arrays(self):
-        """Initialise field arrays on compute device."""
-
-        self.Ex_dev = self.gpuarray.to_gpu(self.Ex)
-        self.Ey_dev = self.gpuarray.to_gpu(self.Ey)
-        self.Ez_dev = self.gpuarray.to_gpu(self.Ez)
-        self.Hx_dev = self.gpuarray.to_gpu(self.Hx)
-        self.Hy_dev = self.gpuarray.to_gpu(self.Hy)
-        self.Hz_dev = self.gpuarray.to_gpu(self.Hz)
-
-    def htod_dispersive_arrays(self):
-        """Initialise dispersive material coefficient arrays on compute device."""
-
-        self.updatecoeffsdispersive_dev = self.gpuarray.to_gpu(self.updatecoeffsdispersive)
-        self.Tx_dev = self.gpuarray.to_gpu(self.Tx)
-        self.Ty_dev = self.gpuarray.to_gpu(self.Ty)
-        self.Tz_dev = self.gpuarray.to_gpu(self.Tz)
-        
-
-class OpenCLGrid(FDTDGrid):
-    """Additional grid methods for solving on compute device using OpenCL."""
-
-    def __init__(self):
-        super().__init__()
-
-        self.clarray = import_module("pyopencl.array")
-
-    def htod_geometry_arrays(self, queue):
-        """Initialise an array for cell edge IDs (ID) on compute device.
-
-        Args:
-            queue: pyopencl queue.
-        """
-
-        self.ID_dev = self.clarray.to_device(queue, self.ID)
-
-    def htod_field_arrays(self, queue):
-        """Initialise field arrays on compute device.
-
-        Args:
-            queue: pyopencl queue.
-        """
-
-        self.Ex_dev = self.clarray.to_device(queue, self.Ex)
-        self.Ey_dev = self.clarray.to_device(queue, self.Ey)
-        self.Ez_dev = self.clarray.to_device(queue, self.Ez)
-        self.Hx_dev = self.clarray.to_device(queue, self.Hx)
-        self.Hy_dev = self.clarray.to_device(queue, self.Hy)
-        self.Hz_dev = self.clarray.to_device(queue, self.Hz)
-
-    def htod_dispersive_arrays(self, queue):
-        """Initialise dispersive material coefficient arrays on compute device.
-
-        Args:
-            queue: pyopencl queue.
-        """
-                 
-        self.updatecoeffsdispersive_dev = self.clarray.to_device(queue, self.updatecoeffsdispersive)
-        # self.updatecoeffsdispersive_dev = self.clarray.to_device(queue, np.ones((95,95,95), dtype=np.float32))
-        self.Tx_dev = self.clarray.to_device(queue, self.Tx)
-        self.Ty_dev = self.clarray.to_device(queue, self.Ty)
-        self.Tz_dev = self.clarray.to_device(queue, self.Tz)
-        
-
 def dispersion_analysis(G):
     """Analysis of numerical dispersion (Taflove et al, 2005, p112) -
         worse case of maximum frequency and minimum wavelength
@@ -444,7 +608,8 @@ def dispersion_analysis(G):
                     try:
                         freqthres = (
                             np.where(
-                                power[freqmaxpower:] < -config.get_model_config().numdispersion["highestfreqthres"]
+                                power[freqmaxpower:]
+                                < -config.get_model_config().numdispersion["highestfreqthres"]
                             )[0][0]
                             + freqmaxpower
                         )
@@ -463,7 +628,8 @@ def dispersion_analysis(G):
                 # If waveform is truncated don't do any further analysis
                 else:
                     results["error"] = (
-                        "waveform does not fit within specified " + "time window and is therefore being truncated."
+                        "waveform does not fit within specified "
+                        + "time window and is therefore being truncated."
                     )
     else:
         results["error"] = "no waveform detected."
@@ -511,7 +677,10 @@ def dispersion_analysis(G):
         results["N"] = minwavelength / delta
 
         # Check grid sampling will result in physical wave propagation
-        if int(np.floor(results["N"])) >= config.get_model_config().numdispersion["mingridsampling"]:
+        if (
+            int(np.floor(results["N"]))
+            >= config.get_model_config().numdispersion["mingridsampling"]
+        ):
             # Numerical phase velocity
             vp = np.pi / (results["N"] * np.arcsin((1 / S) * np.sin((np.pi * S) / results["N"])))
 

gprMax/grid/mpi_grid.py

@@ -0,0 +1,73 @@
+# Copyright (C) 2015-2024: 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/>.
+
+from typing import Optional
+
+from mpi4py import MPI
+
+from gprMax.grid.fdtd_grid import FDTDGrid
+
+
+class MPIGrid(FDTDGrid):
+    def __init__(
+        self,
+        mpi_tasks_x: int,
+        mpi_tasks_y: int,
+        mpi_tasks_z: int,
+        comm: Optional[MPI.Intracomm] = None,
+    ):
+        super().__init__()
+
+        if comm is None:
+            self.comm = MPI.COMM_WORLD
+        else:
+            self.comm = comm
+
+        if mpi_tasks_x * mpi_tasks_y * mpi_tasks_z > self.comm.size:
+            # TODO: Raise expection - insufficient MPI tasks to create the grid as requested
+            pass
+
+        self.mpi_tasks_x = mpi_tasks_x
+        self.mpi_tasks_y = mpi_tasks_y
+        self.mpi_tasks_z = mpi_tasks_z
+
+        self.rank = self.comm.rank
+        self.size = self.comm.size
+
+        self.xmin = 0
+        self.ymin = 0
+        self.zmin = 0
+        self.xmax = 0
+        self.ymax = 0
+        self.zmax = 0
+
+    def initialise_field_arrays(self):
+        super().initialise_field_arrays()
+
+        local_grid_size_x = self.nx // self.mpi_tasks_x
+        local_grid_size_y = self.ny // self.mpi_tasks_y
+        local_grid_size_z = self.nz // self.mpi_tasks_z
+
+        self.xmin = (self.rank % self.nx) * local_grid_size_x
+        self.ymin = ((self.mpi_tasks_x * self.rank) % self.ny) * local_grid_size_y
+        self.zmin = (
+            (self.mpi_tasks_y * self.mpi_tasks_x * self.rank) % self.nz
+        ) * local_grid_size_z
+        self.xmax = self.xmin + local_grid_size_x
+        self.ymax = self.ymin + local_grid_size_y
+        self.zmax = self.zmin + local_grid_size_z

gprMax/grid/opencl_grid.py

@@ -0,0 +1,66 @@
+# Copyright (C) 2015-2024: 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/>.
+
+from importlib import import_module
+
+from gprMax.grid.fdtd_grid import FDTDGrid
+
+
+class OpenCLGrid(FDTDGrid):
+    """Additional grid methods for solving on compute device using OpenCL."""
+
+    def __init__(self):
+        super().__init__()
+
+        self.clarray = import_module("pyopencl.array")
+
+    def htod_geometry_arrays(self, queue):
+        """Initialise an array for cell edge IDs (ID) on compute device.
+
+        Args:
+            queue: pyopencl queue.
+        """
+
+        self.ID_dev = self.clarray.to_device(queue, self.ID)
+
+    def htod_field_arrays(self, queue):
+        """Initialise field arrays on compute device.
+
+        Args:
+            queue: pyopencl queue.
+        """
+
+        self.Ex_dev = self.clarray.to_device(queue, self.Ex)
+        self.Ey_dev = self.clarray.to_device(queue, self.Ey)
+        self.Ez_dev = self.clarray.to_device(queue, self.Ez)
+        self.Hx_dev = self.clarray.to_device(queue, self.Hx)
+        self.Hy_dev = self.clarray.to_device(queue, self.Hy)
+        self.Hz_dev = self.clarray.to_device(queue, self.Hz)
+
+    def htod_dispersive_arrays(self, queue):
+        """Initialise dispersive material coefficient arrays on compute device.
+
+        Args:
+            queue: pyopencl queue.
+        """
+
+        self.updatecoeffsdispersive_dev = self.clarray.to_device(queue, self.updatecoeffsdispersive)
+        # self.updatecoeffsdispersive_dev = self.clarray.to_device(queue, np.ones((95,95,95), dtype=np.float32))
+        self.Tx_dev = self.clarray.to_device(queue, self.Tx)
+        self.Ty_dev = self.clarray.to_device(queue, self.Ty)
+        self.Tz_dev = self.clarray.to_device(queue, self.Tz)

gprMax/model.py

@@ -0,0 +1,237 @@
+# Copyright (C) 2015-2024: 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 datetime
+import itertools
+import logging
+import sys
+
+import humanize
+import numpy as np
+import psutil
+from colorama import Fore, Style, init
+
+from gprMax.grid.cuda_grid import CUDAGrid
+from gprMax.grid.opencl_grid import OpenCLGrid
+
+init()
+
+from terminaltables import SingleTable
+from tqdm import tqdm
+
+import gprMax.config as config
+
+from .cython.yee_cell_build import build_electric_components, build_magnetic_components
+from .fields_outputs import write_hdf5_outputfile
+from .geometry_outputs import save_geometry_views
+from .grid.fdtd_grid import FDTDGrid, dispersion_analysis
+from .materials import process_materials
+from .pml import CFS, build_pml, print_pml_info
+from .snapshots import save_snapshots
+from .utilities.host_info import mem_check_build_all, mem_check_run_all, set_omp_threads
+from .utilities.utilities import get_terminal_width
+
+logger = logging.getLogger(__name__)
+
+
+class Model:
+    """Builds and runs (solves) a model."""
+
+    def __init__(self):
+        self.G = self._create_grid()
+        # Monitor memory usage
+        self.p = None
+
+        # Set number of OpenMP threads to physical threads at this point to be
+        # used with threaded model building methods, e.g. fractals. Can be
+        # changed by the user via #num_threads command in input file or via API
+        # later for use with CPU solver.
+        config.get_model_config().ompthreads = set_omp_threads(config.get_model_config().ompthreads)
+
+    def _create_grid(self) -> FDTDGrid:
+        """Create grid object according to solver.
+
+        Returns:
+            grid: FDTDGrid class describing a grid in a model.
+        """
+        if config.sim_config.general["solver"] == "cpu":
+            grid = FDTDGrid()
+        elif config.sim_config.general["solver"] == "cuda":
+            grid = CUDAGrid()
+        elif config.sim_config.general["solver"] == "opencl":
+            grid = OpenCLGrid()
+
+        return grid
+
+    def build(self):
+        """Builds the Yee cells for a model."""
+
+        G = self.G
+
+        # Monitor memory usage
+        self.p = psutil.Process()
+
+        # Normal model reading/building process; bypassed if geometry information to be reused
+        self.reuse_geometry() if config.get_model_config().reuse_geometry() else self.build_geometry()
+
+        logger.info(
+            f"\nOutput directory: {config.get_model_config().output_file_path.parent.resolve()}"
+        )
+
+        # Adjust position of simple sources and receivers if required
+        model_num = config.sim_config.current_model
+        G.update_simple_source_positions(step=model_num)
+        G.update_receiver_positions(step=model_num)
+
+        # Write files for any geometry views and geometry object outputs
+        gvs = G.geometryviews + [gv for sg in G.subgrids for gv in sg.geometryviews]
+        if not gvs and not G.geometryobjectswrite and config.sim_config.args.geometry_only:
+            logger.exception("\nNo geometry views or geometry objects found.")
+            raise ValueError
+        save_geometry_views(gvs)
+
+        if G.geometryobjectswrite:
+            logger.info("")
+            for i, go in enumerate(G.geometryobjectswrite):
+                pbar = tqdm(
+                    total=go.datawritesize,
+                    unit="byte",
+                    unit_scale=True,
+                    desc=f"Writing geometry object file {i + 1}/{len(G.geometryobjectswrite)}, "
+                    + f"{go.filename_hdf5.name}",
+                    ncols=get_terminal_width() - 1,
+                    file=sys.stdout,
+                    disable=not config.sim_config.general["progressbars"],
+                )
+                go.write_hdf5(G, pbar)
+                pbar.close()
+            logger.info("")
+
+    def build_geometry(self):
+        logger.info(config.get_model_config().inputfilestr)
+        self.G.build()
+
+    def reuse_geometry(self):
+        s = (
+            f"\n--- Model {config.get_model_config().appendmodelnumber}/{config.sim_config.model_end}, "
+            f"input file (not re-processed, i.e. geometry fixed): "
+            f"{config.sim_config.input_file_path}"
+        )
+        config.get_model_config().inputfilestr = (
+            Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
+        )
+        logger.basic(config.get_model_config().inputfilestr)
+        for grid in [self.G] + self.G.subgrids:
+            grid.iteration = 0  # Reset current iteration number
+            grid.reset_fields()
+
+    def write_output_data(self):
+        """Writes output data, i.e. field data for receivers and snapshots to
+        file(s).
+        """
+
+        # Write output data to file if they are any receivers in any grids
+        sg_rxs = [True for sg in self.G.subgrids if sg.rxs]
+        sg_tls = [True for sg in self.G.subgrids if sg.transmissionlines]
+        if self.G.rxs or sg_rxs or self.G.transmissionlines or sg_tls:
+            write_hdf5_outputfile(config.get_model_config().output_file_path_ext, self.G)
+
+        # Write any snapshots to file for each grid
+        for grid in [self.G] + self.G.subgrids:
+            if grid.snapshots:
+                save_snapshots(grid)
+
+    def solve(self, solver):
+        """Solve using FDTD method.
+
+        Args:
+            solver: solver object.
+        """
+
+        # Print information about and check OpenMP threads
+        if config.sim_config.general["solver"] == "cpu":
+            logger.basic(
+                f"\nModel {config.sim_config.current_model + 1}/{config.sim_config.model_end} "
+                f"on {config.sim_config.hostinfo['hostname']} "
+                f"with OpenMP backend using {config.get_model_config().ompthreads} thread(s)"
+            )
+            if config.get_model_config().ompthreads > config.sim_config.hostinfo["physicalcores"]:
+                logger.warning(
+                    f"You have specified more threads ({config.get_model_config().ompthreads}) "
+                    f"than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). "
+                    f"This may lead to degraded performance."
+                )
+        elif config.sim_config.general["solver"] in ["cuda", "opencl"]:
+            if config.sim_config.general["solver"] == "opencl":
+                solvername = "OpenCL"
+                platformname = (
+                    " ".join(config.get_model_config().device["dev"].platform.name.split())
+                    + " with "
+                )
+                devicename = (
+                    f'Device {config.get_model_config().device["deviceID"]}: '
+                    f'{" ".join(config.get_model_config().device["dev"].name.split())}'
+                )
+            else:
+                solvername = "CUDA"
+                platformname = ""
+                devicename = (
+                    f'Device {config.get_model_config().device["deviceID"]}: '
+                    f'{" ".join(config.get_model_config().device["dev"].name().split())}'
+                )
+
+            logger.basic(
+                f"\nModel {config.sim_config.current_model + 1}/{config.sim_config.model_end} "
+                f"solving on {config.sim_config.hostinfo['hostname']} "
+                f"with {solvername} backend using {platformname}{devicename}"
+            )
+
+        # Prepare iterator
+        if config.sim_config.general["progressbars"]:
+            iterator = tqdm(
+                range(self.G.iterations),
+                desc="|--->",
+                ncols=get_terminal_width() - 1,
+                file=sys.stdout,
+                disable=not config.sim_config.general["progressbars"],
+            )
+        else:
+            iterator = range(self.G.iterations)
+
+        # Run solver
+        solver.solve(iterator)
+
+        # Write output data, i.e. field data for receivers and snapshots to file(s)
+        self.write_output_data()
+
+        # Print information about memory usage and solving time for a model
+        # Add a string on device (GPU) memory usage if applicable
+        mem_str = ""
+        if config.sim_config.general["solver"] == "cuda":
+            mem_str = f" host + ~{humanize.naturalsize(solver.memused)} device"
+        elif config.sim_config.general["solver"] == "opencl":
+            mem_str = f" host + unknown for device"
+
+        logger.info(
+            f"\nMemory used (estimated): "
+            + f"~{humanize.naturalsize(self.p.memory_full_info().uss)}{mem_str}"
+        )
+        logger.info(
+            f"Time taken: "
+            + f"{humanize.precisedelta(datetime.timedelta(seconds=solver.solvetime), format='%0.4f')}"
+        )

gprMax/model_build_run.py

@@ -1,420 +0,0 @@
-# Copyright (C) 2015-2024: 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 datetime
-import itertools
-import logging
-import sys
-
-import humanize
-import numpy as np
-import psutil
-from colorama import Fore, Style, init
-
-init()
-
-from terminaltables import SingleTable
-from tqdm import tqdm
-
-import gprMax.config as config
-
-from .cython.yee_cell_build import build_electric_components, build_magnetic_components
-from .fields_outputs import write_hdf5_outputfile
-from .geometry_outputs import save_geometry_views
-from .grid import dispersion_analysis
-from .hash_cmds_file import parse_hash_commands
-from .materials import process_materials
-from .pml import CFS, build_pml, print_pml_info
-from .scene import Scene
-from .snapshots import save_snapshots
-from .utilities.host_info import mem_check_build_all, mem_check_run_all, set_omp_threads
-from .utilities.utilities import get_terminal_width
-
-logger = logging.getLogger(__name__)
-
-
-class ModelBuildRun:
-    """Builds and runs (solves) a model."""
-
-    def __init__(self, G):
-        self.G = G
-        # Monitor memory usage
-        self.p = None
-
-        # Set number of OpenMP threads to physical threads at this point to be
-        # used with threaded model building methods, e.g. fractals. Can be
-        # changed by the user via #num_threads command in input file or via API
-        # later for use with CPU solver.
-        config.get_model_config().ompthreads = set_omp_threads(config.get_model_config().ompthreads)
-
-    def build(self):
-        """Builds the Yee cells for a model."""
-
-        G = self.G
-
-        # Monitor memory usage
-        self.p = psutil.Process()
-
-        # Normal model reading/building process; bypassed if geometry information to be reused
-        self.reuse_geometry() if config.get_model_config().reuse_geometry else self.build_geometry()
-
-        logger.info(f"\nOutput directory: {config.get_model_config().output_file_path.parent.resolve()}")
-
-        # 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 config.model_num == 0:
-                    if (
-                        source.xcoord + G.srcsteps[0] * config.sim_config.model_end < 0
-                        or source.xcoord + G.srcsteps[0] * config.sim_config.model_end > G.nx
-                        or source.ycoord + G.srcsteps[1] * config.sim_config.model_end < 0
-                        or source.ycoord + G.srcsteps[1] * config.sim_config.model_end > G.ny
-                        or source.zcoord + G.srcsteps[2] * config.sim_config.model_end < 0
-                        or source.zcoord + G.srcsteps[2] * config.sim_config.model_end > G.nz
-                    ):
-                        logger.exception("Source(s) will be stepped to a position outside the domain.")
-                        raise ValueError
-                source.xcoord = source.xcoordorigin + config.model_num * G.srcsteps[0]
-                source.ycoord = source.ycoordorigin + config.model_num * G.srcsteps[1]
-                source.zcoord = source.zcoordorigin + config.model_num * G.srcsteps[2]
-        if G.rxsteps[0] != 0 or G.rxsteps[1] != 0 or G.rxsteps[2] != 0:
-            for receiver in G.rxs:
-                if config.model_num == 0:
-                    if (
-                        receiver.xcoord + G.rxsteps[0] * config.sim_config.model_end < 0
-                        or receiver.xcoord + G.rxsteps[0] * config.sim_config.model_end > G.nx
-                        or receiver.ycoord + G.rxsteps[1] * config.sim_config.model_end < 0
-                        or receiver.ycoord + G.rxsteps[1] * config.sim_config.model_end > G.ny
-                        or receiver.zcoord + G.rxsteps[2] * config.sim_config.model_end < 0
-                        or receiver.zcoord + G.rxsteps[2] * config.sim_config.model_end > G.nz
-                    ):
-                        logger.exception("Receiver(s) will be stepped to a position outside the domain.")
-                        raise ValueError
-                receiver.xcoord = receiver.xcoordorigin + config.model_num * G.rxsteps[0]
-                receiver.ycoord = receiver.ycoordorigin + config.model_num * G.rxsteps[1]
-                receiver.zcoord = receiver.zcoordorigin + config.model_num * G.rxsteps[2]
-
-        # Write files for any geometry views and geometry object outputs
-        gvs = G.geometryviews + [gv for sg in G.subgrids for gv in sg.geometryviews]
-        if not gvs and not G.geometryobjectswrite and config.sim_config.args.geometry_only:
-            logger.exception("\nNo geometry views or geometry objects found.")
-            raise ValueError
-        save_geometry_views(gvs)
-
-        if G.geometryobjectswrite:
-            logger.info("")
-            for i, go in enumerate(G.geometryobjectswrite):
-                pbar = tqdm(
-                    total=go.datawritesize,
-                    unit="byte",
-                    unit_scale=True,
-                    desc=f"Writing geometry object file {i + 1}/{len(G.geometryobjectswrite)}, "
-                    + f"{go.filename_hdf5.name}",
-                    ncols=get_terminal_width() - 1,
-                    file=sys.stdout,
-                    disable=not config.sim_config.general["progressbars"],
-                )
-                go.write_hdf5(G, pbar)
-                pbar.close()
-            logger.info("")
-
-    def build_geometry(self):
-        G = self.G
-
-        logger.info(config.get_model_config().inputfilestr)
-
-        # Build objects in the scene and check memory for building
-        self.build_scene()
-
-        # Print info on any subgrids
-        for sg in G.subgrids:
-            sg.print_info()
-
-        # Combine available grids
-        grids = [G] + G.subgrids
-
-        # Check for dispersive materials (and specific type)
-        for grid in grids:
-            if config.get_model_config().materials["maxpoles"] != 0:
-                config.get_model_config().materials["drudelorentz"] = any(
-                    [m for m in grid.materials if "drude" in m.type or "lorentz" in m.type]
-                )
-
-        # Set data type if any dispersive materials (must be done before memory checks)
-        if config.get_model_config().materials["maxpoles"] != 0:
-            config.get_model_config().set_dispersive_material_types()
-
-        # Check memory requirements to build model/scene (different to memory
-        # requirements to run model when FractalVolumes/FractalSurfaces are
-        # used as these can require significant additional memory)
-        total_mem_build, mem_strs_build = mem_check_build_all(grids)
-
-        # Check memory requirements to run model
-        total_mem_run, mem_strs_run = mem_check_run_all(grids)
-
-        if total_mem_build > total_mem_run:
-            logger.info(
-                f'\nMemory required (estimated): {" + ".join(mem_strs_build)} + '
-                f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
-                f"overhead = {humanize.naturalsize(total_mem_build)}"
-            )
-        else:
-            logger.info(
-                f'\nMemory required (estimated): {" + ".join(mem_strs_run)} + '
-                f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
-                f"overhead = {humanize.naturalsize(total_mem_run)}"
-            )
-
-        # Build grids
-        gridbuilders = [GridBuilder(grid) for grid in grids]
-        for gb in gridbuilders:
-            # Set default CFS parameter for PMLs if not user provided
-            if not gb.grid.pmls["cfs"]:
-                gb.grid.pmls["cfs"] = [CFS()]
-            logger.info(print_pml_info(gb.grid))
-            if not all(value == 0 for value in gb.grid.pmls["thickness"].values()):
-                gb.build_pmls()
-            if gb.grid.averagevolumeobjects:
-                gb.build_components()
-            gb.tm_grid_update()
-            gb.update_voltage_source_materials()
-            gb.grid.initialise_field_arrays()
-            gb.grid.initialise_std_update_coeff_arrays()
-            if config.get_model_config().materials["maxpoles"] > 0:
-                gb.grid.initialise_dispersive_arrays()
-                gb.grid.initialise_dispersive_update_coeff_array()
-            gb.build_materials()
-
-            # Check to see if numerical dispersion might be a problem
-            results = dispersion_analysis(gb.grid)
-            if results["error"]:
-                logger.warning(
-                    f"\nNumerical dispersion analysis [{gb.grid.name}] " f"not carried out as {results['error']}"
-                )
-            elif results["N"] < config.get_model_config().numdispersion["mingridsampling"]:
-                logger.exception(
-                    f"\nNon-physical wave propagation in [{gb.grid.name}] "
-                    f"detected. Material '{results['material'].ID}' "
-                    f"has wavelength sampled by {results['N']} cells, "
-                    f"less than required minimum for physical wave "
-                    f"propagation. Maximum significant frequency "
-                    f"estimated as {results['maxfreq']:g}Hz"
-                )
-                raise ValueError
-            elif (
-                results["deltavp"]
-                and np.abs(results["deltavp"]) > config.get_model_config().numdispersion["maxnumericaldisp"]
-            ):
-                logger.warning(
-                    f"\n[{gb.grid.name}] has potentially significant "
-                    f"numerical dispersion. Estimated largest physical "
-                    f"phase-velocity error is {results['deltavp']:.2f}% "
-                    f"in material '{results['material'].ID}' whose "
-                    f"wavelength sampled by {results['N']} cells. "
-                    f"Maximum significant frequency estimated as "
-                    f"{results['maxfreq']:g}Hz"
-                )
-            elif results["deltavp"]:
-                logger.info(
-                    f"\nNumerical dispersion analysis [{gb.grid.name}]: "
-                    f"estimated largest physical phase-velocity error is "
-                    f"{results['deltavp']:.2f}% in material '{results['material'].ID}' "
-                    f"whose wavelength sampled by {results['N']} cells. "
-                    f"Maximum significant frequency estimated as "
-                    f"{results['maxfreq']:g}Hz"
-                )
-
-    def reuse_geometry(self):
-        s = (
-            f"\n--- Model {config.get_model_config().appendmodelnumber}/{config.sim_config.model_end}, "
-            f"input file (not re-processed, i.e. geometry fixed): "
-            f"{config.sim_config.input_file_path}"
-        )
-        config.get_model_config().inputfilestr = (
-            Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
-        )
-        logger.basic(config.get_model_config().inputfilestr)
-        for grid in [self.G] + self.G.subgrids:
-            grid.iteration = 0  # Reset current iteration number
-            grid.reset_fields()
-
-    def build_scene(self):
-        # API for multiple scenes / model runs
-        scene = config.get_model_config().get_scene()
-
-        # If there is no scene, process the hash commands
-        if not scene:
-            scene = Scene()
-            config.sim_config.scenes.append(scene)
-            # Parse the input file into user objects and add them to the scene
-            scene = parse_hash_commands(scene)
-
-        # Creates the internal simulation objects
-        scene.create_internal_objects(self.G)
-
-        return scene
-
-    def write_output_data(self):
-        """Writes output data, i.e. field data for receivers and snapshots to
-        file(s).
-        """
-
-        # Write output data to file if they are any receivers in any grids
-        sg_rxs = [True for sg in self.G.subgrids if sg.rxs]
-        sg_tls = [True for sg in self.G.subgrids if sg.transmissionlines]
-        if self.G.rxs or sg_rxs or self.G.transmissionlines or sg_tls:
-            write_hdf5_outputfile(config.get_model_config().output_file_path_ext, self.G)
-
-        # Write any snapshots to file for each grid
-        for grid in [self.G] + self.G.subgrids:
-            if grid.snapshots:
-                save_snapshots(grid)
-
-    def solve(self, solver):
-        """Solve using FDTD method.
-
-        Args:
-            solver: solver object.
-        """
-
-        # Print information about and check OpenMP threads
-        if config.sim_config.general["solver"] == "cpu":
-            logger.basic(
-                f"\nModel {config.model_num + 1}/{config.sim_config.model_end} "
-                f"on {config.sim_config.hostinfo['hostname']} "
-                f"with OpenMP backend using {config.get_model_config().ompthreads} thread(s)"
-            )
-            if config.get_model_config().ompthreads > config.sim_config.hostinfo["physicalcores"]:
-                logger.warning(
-                    f"You have specified more threads ({config.get_model_config().ompthreads}) "
-                    f"than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). "
-                    f"This may lead to degraded performance."
-                )
-        elif config.sim_config.general["solver"] in ["cuda", "opencl"]:
-            if config.sim_config.general["solver"] == "opencl":
-                solvername = "OpenCL"
-                platformname = " ".join(config.get_model_config().device["dev"].platform.name.split()) + " with "
-                devicename = (
-                    f'Device {config.get_model_config().device["deviceID"]}: '
-                    f'{" ".join(config.get_model_config().device["dev"].name.split())}'
-                )
-            else:
-                solvername = "CUDA"
-                platformname = ""
-                devicename = (
-                    f'Device {config.get_model_config().device["deviceID"]}: '
-                    f'{" ".join(config.get_model_config().device["dev"].name().split())}'
-                )
-
-            logger.basic(
-                f"\nModel {config.model_num + 1}/{config.sim_config.model_end} "
-                f"solving on {config.sim_config.hostinfo['hostname']} "
-                f"with {solvername} backend using {platformname}{devicename}"
-            )
-
-        # Prepare iterator
-        if config.sim_config.general["progressbars"]:
-            iterator = tqdm(
-                range(self.G.iterations),
-                desc="|--->",
-                ncols=get_terminal_width() - 1,
-                file=sys.stdout,
-                disable=not config.sim_config.general["progressbars"],
-            )
-        else:
-            iterator = range(self.G.iterations)
-
-        # Run solver
-        solver.solve(iterator)
-
-        # Write output data, i.e. field data for receivers and snapshots to file(s)
-        self.write_output_data()
-
-        # Print information about memory usage and solving time for a model
-        # Add a string on device (GPU) memory usage if applicable
-        mem_str = ""
-        if config.sim_config.general["solver"] == "cuda":
-            mem_str = f" host + ~{humanize.naturalsize(solver.memused)} device"
-        elif config.sim_config.general["solver"] == "opencl":
-            mem_str = f" host + unknown for device"      
-
-        logger.info(f"\nMemory used (estimated): " + f"~{humanize.naturalsize(self.p.memory_full_info().uss)}{mem_str}")
-        logger.info(
-            f"Time taken: " + f"{humanize.precisedelta(datetime.timedelta(seconds=solver.solvetime), format='%0.4f')}"
-        )
-
-
-class GridBuilder:
-    def __init__(self, grid):
-        self.grid = grid
-
-    def build_pmls(self):
-        pbar = tqdm(
-            total=sum(1 for value in self.grid.pmls["thickness"].values() if value > 0),
-            desc=f"Building PML boundaries [{self.grid.name}]",
-            ncols=get_terminal_width() - 1,
-            file=sys.stdout,
-            disable=not config.sim_config.general["progressbars"],
-        )
-        for pml_id, thickness in self.grid.pmls["thickness"].items():
-            if thickness > 0:
-                build_pml(self.grid, pml_id, thickness)
-                pbar.update()
-        pbar.close()
-
-    def build_components(self):
-        # Build the model, i.e. set the material properties (ID) for every edge
-        # of every Yee cell
-        logger.info("")
-        pbar = tqdm(
-            total=2,
-            desc=f"Building Yee cells [{self.grid.name}]",
-            ncols=get_terminal_width() - 1,
-            file=sys.stdout,
-            disable=not config.sim_config.general["progressbars"],
-        )
-        build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid)
-        pbar.update()
-        build_magnetic_components(self.grid.solid, self.grid.rigidH, self.grid.ID, self.grid)
-        pbar.update()
-        pbar.close()
-
-    def tm_grid_update(self):
-        if config.get_model_config().mode == "2D TMx":
-            self.grid.tmx()
-        elif config.get_model_config().mode == "2D TMy":
-            self.grid.tmy()
-        elif config.get_model_config().mode == "2D TMz":
-            self.grid.tmz()
-
-    def update_voltage_source_materials(self):
-        # Process any voltage sources (that have resistance) to create a new
-        # material at the source location
-        for voltagesource in self.grid.voltagesources:
-            voltagesource.create_material(self.grid)
-
-    def build_materials(self):
-        # Process complete list of materials - calculate update coefficients,
-        # store in arrays, and build text list of materials/properties
-        materialsdata = process_materials(self.grid)
-        materialstable = SingleTable(materialsdata)
-        materialstable.outer_border = False
-        materialstable.justify_columns[0] = "right"
-
-        logger.info(f"\nMaterials [{self.grid.name}]:")
-        logger.info(materialstable.table)

gprMax/solvers.py

@@ -18,74 +18,21 @@
 
 import gprMax.config as config
 
-from .grid import CUDAGrid, FDTDGrid, OpenCLGrid
+from .grid.cuda_grid import CUDAGrid
+from .grid.fdtd_grid import FDTDGrid
+from .grid.opencl_grid import OpenCLGrid
+from .subgrids.updates import SubgridUpdates
 from .subgrids.updates import create_updates as create_subgrid_updates
-from .updates.updates import CPUUpdates, CUDAUpdates, OpenCLUpdates
+from .updates.cpu_updates import CPUUpdates
-
+from .updates.cuda_updates import CUDAUpdates
-
+from .updates.opencl_updates import OpenCLUpdates
-def create_G():
+from .updates.updates import Updates
-    """Create grid object according to solver.
-
-    Returns:
-        G: FDTDGrid class describing a grid in a model.
-    """
-
-    if config.sim_config.general["solver"] == "cpu":
-        G = FDTDGrid()
-    elif config.sim_config.general["solver"] == "cuda":
-        G = CUDAGrid()
-    elif config.sim_config.general["solver"] == "opencl":
-        G = OpenCLGrid()
-
-    return G
-
-
-def create_solver(G):
-    """Create configured solver object.
-
-    N.B. A large range of different functions exist to advance the time step for
-            dispersive materials. The correct function is set by the
-            set_dispersive_updates method, based on the required numerical
-            precision and dispersive material type.
-            This is done for solvers running on CPU, i.e. where Cython is used.
-            CUDA and OpenCL dispersive material functions are handled through
-            templating and substitution at runtime.
-
-    Args:
-        G: FDTDGrid class describing a grid in a model.
-
-    Returns:
-        solver: Solver object.
-    """
-
-    if config.sim_config.general["subgrid"]:
-        updates = create_subgrid_updates(G)
-        if config.get_model_config().materials["maxpoles"] != 0:
-            # Set dispersive update functions for both SubgridUpdates and
-            # SubgridUpdaters subclasses
-            updates.set_dispersive_updates()
-            for u in updates.updaters:
-                u.set_dispersive_updates()
-        solver = Solver(updates, hsg=True)
-    elif config.sim_config.general["solver"] == "cpu":
-        updates = CPUUpdates(G)
-        if config.get_model_config().materials["maxpoles"] != 0:
-            updates.set_dispersive_updates()
-        solver = Solver(updates)
-    elif config.sim_config.general["solver"] == "cuda":
-        updates = CUDAUpdates(G)
-        solver = Solver(updates)
-    elif config.sim_config.general["solver"] == "opencl":
-        updates = OpenCLUpdates(G)
-        solver = Solver(updates)
-
-    return solver
 
 
 class Solver:
     """Generic solver for Update objects"""
 
-    def __init__(self, updates, hsg=False):
+    def __init__(self, updates: Updates):
         """
         Args:
             updates: Updates contains methods to run FDTD algorithm.
@@ -93,7 +40,6 @@ class Solver:
         """
 
         self.updates = updates
-        self.hsg = hsg
         self.solvetime = 0
         self.memused = 0
 
@@ -112,17 +58,57 @@ class Solver:
             self.updates.update_magnetic()
             self.updates.update_magnetic_pml()
             self.updates.update_magnetic_sources()
-            if self.hsg:
+            if isinstance(self.updates, SubgridUpdates):
                 self.updates.hsg_2()
             self.updates.update_electric_a()
             self.updates.update_electric_pml()
             self.updates.update_electric_sources()
-            if self.hsg:
+            if isinstance(self.updates, SubgridUpdates):
                 self.updates.hsg_1()
             self.updates.update_electric_b()
-            if config.sim_config.general["solver"] == "cuda":
+            if isinstance(self.updates, CUDAUpdates):
                 self.memused = self.updates.calculate_memory_used(iteration)
 
         self.updates.finalise()
         self.solvetime = self.updates.calculate_solve_time()
         self.updates.cleanup()
+
+
+def create_solver(grid: FDTDGrid) -> Solver:
+    """Create configured solver object.
+
+    N.B. A large range of different functions exist to advance the time
+    step for dispersive materials. The correct function is set by the
+    set_dispersive_updates method, based on the required numerical
+    precision and dispersive material type. This is done for solvers
+    running on CPU, i.e. where Cython is used. CUDA and OpenCL
+    dispersive material functions are handled through templating and
+    substitution at runtime.
+
+    Args:
+        G: FDTDGrid class describing a grid in a model.
+
+    Returns:
+        solver: Solver object.
+    """
+
+    if config.sim_config.general["subgrid"]:
+        updates = create_subgrid_updates(grid)
+        if config.get_model_config().materials["maxpoles"] != 0:
+            # Set dispersive update functions for both SubgridUpdates and
+            # SubgridUpdaters subclasses
+            updates.set_dispersive_updates()
+            for u in updates.updaters:
+                u.set_dispersive_updates()
+    elif type(grid) is FDTDGrid:
+        updates = CPUUpdates(grid)
+        if config.get_model_config().materials["maxpoles"] != 0:
+            updates.set_dispersive_updates()
+    elif type(grid) is CUDAGrid:
+        updates = CUDAUpdates(grid)
+    elif type(grid) is OpenCLGrid:
+        updates = OpenCLUpdates(grid)
+
+    solver = Solver(updates)
+
+    return solver

gprMax/subgrids/grid.py

@@ -17,8 +17,9 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import logging
+from abc import abstractmethod
 
-from ..grid import FDTDGrid
+from gprMax.grid.fdtd_grid import FDTDGrid
 
 logger = logging.getLogger(__name__)
 
@@ -62,3 +63,23 @@ class SubGridBaseGrid(FDTDGrid):
         self.n_boundary_cells_z = d_to_pml + self.pmls["thickness"]["z0"]
 
         self.interpolation = kwargs["interpolation"]
+
+    @abstractmethod
+    def update_magnetic_is(self, precursors):
+        pass
+
+    @abstractmethod
+    def update_electric_is(self, precursors):
+        pass
+
+    @abstractmethod
+    def update_electric_os(self, main_grid):
+        pass
+
+    @abstractmethod
+    def update_magnetic_os(self, main_grid):
+        pass
+
+    @abstractmethod
+    def print_info(self):
+        pass

gprMax/subgrids/updates.py

@@ -18,7 +18,7 @@
 
 import logging
 
-from ..updates.updates import CPUUpdates
+from ..updates.cpu_updates import CPUUpdates
 from .precursor_nodes import PrecursorNodes, PrecursorNodesFiltered
 from .subgrid_hsg import SubGridHSG
 

gprMax/taskfarm.py

@@ -51,7 +51,7 @@ class TaskfarmExecutor(object):
     `gprMax` models in parallel is given below.
     >>> from mpi4py import MPI
     >>> from gprMax.taskfarm import TaskfarmExecutor
-    >>> from gprMax.model_build_run import run_model
+    >>> from gprMax.model import run_model
     >>> # choose an MPI.Intracomm for communication (MPI.COMM_WORLD by default)
     >>> comm = MPI.COMM_WORLD
     >>> # choose a target function
@@ -269,7 +269,9 @@ class TaskfarmExecutor(object):
             for i, worker in enumerate(self.workers):
                 if self.comm.Iprobe(source=worker, tag=Tags.DONE):
                     job_idx, result = self.comm.recv(source=worker, tag=Tags.DONE)
-                    logger.debug(f"({self.comm.name}) - Received finished job {job_idx} from worker {worker:d}.")
+                    logger.debug(
+                        f"({self.comm.name}) - Received finished job {job_idx} from worker {worker:d}."
+                    )
                     results[job_idx] = result
                     self.busy[i] = False
                 elif self.comm.Iprobe(source=worker, tag=Tags.READY):
@@ -277,7 +279,9 @@ class TaskfarmExecutor(object):
                         self.comm.recv(source=worker, tag=Tags.READY)
                         self.busy[i] = True
                         job_idx = num_jobs - len(my_jobs)
-                        logger.debug(f"({self.comm.name}) - Sending job {job_idx} to worker {worker:d}.")
+                        logger.debug(
+                            f"({self.comm.name}) - Sending job {job_idx} to worker {worker:d}."
+                        )
                         self.comm.send((job_idx, my_jobs.pop(0)), dest=worker, tag=Tags.START)
                 elif self.comm.Iprobe(source=worker, tag=Tags.EXIT):
                     logger.debug(f"({self.comm.name}) - Worker on rank {worker:d} has terminated.")

gprMax/updates/cpu_updates.py

@@ -0,0 +1,213 @@
+# Copyright (C) 2015-2024: 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/>.
+
+from importlib import import_module
+
+from gprMax import config
+from gprMax.cython.fields_updates_normal import update_electric as update_electric_cpu
+from gprMax.cython.fields_updates_normal import update_magnetic as update_magnetic_cpu
+from gprMax.fields_outputs import store_outputs as store_outputs_cpu
+from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.updates.updates import Updates
+from gprMax.utilities.utilities import timer
+
+
+class CPUUpdates(Updates):
+    """Defines update functions for CPU-based solver."""
+
+    def __init__(self, G: FDTDGrid):
+        """
+        Args:
+            G: FDTDGrid class describing a grid in a model.
+        """
+
+        self.grid = G
+
+    def store_outputs(self):
+        """Stores field component values for every receiver and transmission line."""
+        store_outputs_cpu(self.grid)
+
+    def store_snapshots(self, iteration):
+        """Stores any snapshots.
+
+        Args:
+            iteration: int for iteration number.
+        """
+        for snap in self.grid.snapshots:
+            if snap.time == iteration + 1:
+                snap.store(self.grid)
+
+    def update_magnetic(self):
+        """Updates magnetic field components."""
+        update_magnetic_cpu(
+            self.grid.nx,
+            self.grid.ny,
+            self.grid.nz,
+            config.get_model_config().ompthreads,
+            self.grid.updatecoeffsH,
+            self.grid.ID,
+            self.grid.Ex,
+            self.grid.Ey,
+            self.grid.Ez,
+            self.grid.Hx,
+            self.grid.Hy,
+            self.grid.Hz,
+        )
+
+    def update_magnetic_pml(self):
+        """Updates magnetic field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_magnetic()
+
+    def update_magnetic_sources(self):
+        """Updates magnetic field components from sources."""
+        for source in self.grid.transmissionlines + self.grid.magneticdipoles:
+            source.update_magnetic(
+                self.grid.iteration,
+                self.grid.updatecoeffsH,
+                self.grid.ID,
+                self.grid.Hx,
+                self.grid.Hy,
+                self.grid.Hz,
+                self.grid,
+            )
+
+    def update_electric_a(self):
+        """Updates electric field components."""
+        # All materials are non-dispersive so do standard update.
+        if config.get_model_config().materials["maxpoles"] == 0:
+            update_electric_cpu(
+                self.grid.nx,
+                self.grid.ny,
+                self.grid.nz,
+                config.get_model_config().ompthreads,
+                self.grid.updatecoeffsE,
+                self.grid.ID,
+                self.grid.Ex,
+                self.grid.Ey,
+                self.grid.Ez,
+                self.grid.Hx,
+                self.grid.Hy,
+                self.grid.Hz,
+            )
+
+        # 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).
+        else:
+            self.dispersive_update_a(
+                self.grid.nx,
+                self.grid.ny,
+                self.grid.nz,
+                config.get_model_config().ompthreads,
+                config.get_model_config().materials["maxpoles"],
+                self.grid.updatecoeffsE,
+                self.grid.updatecoeffsdispersive,
+                self.grid.ID,
+                self.grid.Tx,
+                self.grid.Ty,
+                self.grid.Tz,
+                self.grid.Ex,
+                self.grid.Ey,
+                self.grid.Ez,
+                self.grid.Hx,
+                self.grid.Hy,
+                self.grid.Hz,
+            )
+
+    def update_electric_pml(self):
+        """Updates electric field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_electric()
+
+    def update_electric_sources(self):
+        """Updates electric field components from sources -
+        update any Hertzian dipole sources last.
+        """
+        for source in (
+            self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles
+        ):
+            source.update_electric(
+                self.grid.iteration,
+                self.grid.updatecoeffsE,
+                self.grid.ID,
+                self.grid.Ex,
+                self.grid.Ey,
+                self.grid.Ez,
+                self.grid,
+            )
+        self.grid.iteration += 1
+
+    def update_electric_b(self):
+        """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 config.get_model_config().materials["maxpoles"] > 0:
+            self.dispersive_update_b(
+                self.grid.nx,
+                self.grid.ny,
+                self.grid.nz,
+                config.get_model_config().ompthreads,
+                config.get_model_config().materials["maxpoles"],
+                self.grid.updatecoeffsdispersive,
+                self.grid.ID,
+                self.grid.Tx,
+                self.grid.Ty,
+                self.grid.Tz,
+                self.grid.Ex,
+                self.grid.Ey,
+                self.grid.Ez,
+            )
+
+    def set_dispersive_updates(self):
+        """Sets dispersive update functions."""
+
+        poles = "multi" if config.get_model_config().materials["maxpoles"] > 1 else "1"
+        precision = "float" if config.sim_config.general["precision"] == "single" else "double"
+        dispersion = (
+            "complex"
+            if config.get_model_config().materials["dispersivedtype"]
+            == config.sim_config.dtypes["complex"]
+            else "real"
+        )
+
+        update_f = "update_electric_dispersive_{}pole_{}_{}_{}"
+        disp_a = update_f.format(poles, "A", precision, dispersion)
+        disp_b = update_f.format(poles, "B", precision, dispersion)
+
+        disp_a_f = getattr(import_module("gprMax.cython.fields_updates_dispersive"), disp_a)
+        disp_b_f = getattr(import_module("gprMax.cython.fields_updates_dispersive"), disp_b)
+
+        self.dispersive_update_a = disp_a_f
+        self.dispersive_update_b = disp_b_f
+
+    def time_start(self):
+        """Starts timer used to calculate solving time for model."""
+        self.timestart = timer()
+
+    def calculate_solve_time(self):
+        """Calculates solving time for model."""
+        return timer() - self.timestart
+
+    def finalise(self):
+        pass
+
+    def cleanup(self):
+        pass

gprMax/updates/cuda_updates.py

@@ -0,0 +1,629 @@
+# Copyright (C) 2015-2024: 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 importlib import import_module
+
+import humanize
+import numpy as np
+from jinja2 import Environment, PackageLoader
+
+from gprMax import config
+from gprMax.cuda_opencl import (
+    knl_fields_updates,
+    knl_snapshots,
+    knl_source_updates,
+    knl_store_outputs,
+)
+from gprMax.grid.cuda_grid import CUDAGrid
+from gprMax.receivers import dtoh_rx_array, htod_rx_arrays
+from gprMax.snapshots import Snapshot, dtoh_snapshot_array, htod_snapshot_array
+from gprMax.sources import htod_src_arrays
+from gprMax.updates.updates import Updates
+from gprMax.utilities.utilities import round32
+
+logger = logging.getLogger(__name__)
+
+
+class CUDAUpdates(Updates):
+    """Defines update functions for GPU-based (CUDA) solver."""
+
+    def __init__(self, G: CUDAGrid):
+        """
+        Args:
+            G: CUDAGrid class describing a grid in a model.
+        """
+
+        self.grid = G
+
+        # Import PyCUDA modules
+        self.drv = import_module("pycuda.driver")
+        self.source_module = getattr(import_module("pycuda.compiler"), "SourceModule")
+        self.drv.init()
+
+        # Create device handle and context on specific GPU device (and make it current context)
+        self.dev = config.get_model_config().device["dev"]
+        self.ctx = self.dev.make_context()
+
+        # Set common substitutions for use in kernels
+        # Substitutions in function arguments
+        self.subs_name_args = {
+            "REAL": config.sim_config.dtypes["C_float_or_double"],
+            "COMPLEX": config.get_model_config().materials["dispersiveCdtype"],
+        }
+        # Substitutions in function bodies
+        self.subs_func = {
+            "REAL": config.sim_config.dtypes["C_float_or_double"],
+            "CUDA_IDX": "int i = blockIdx.x * blockDim.x + threadIdx.x;",
+            "NX_FIELDS": self.grid.nx + 1,
+            "NY_FIELDS": self.grid.ny + 1,
+            "NZ_FIELDS": self.grid.nz + 1,
+            "NX_ID": self.grid.ID.shape[1],
+            "NY_ID": self.grid.ID.shape[2],
+            "NZ_ID": self.grid.ID.shape[3],
+        }
+
+        # Enviroment for templating kernels
+        self.env = Environment(loader=PackageLoader("gprMax", "cuda_opencl"))
+
+        # Initialise arrays on GPU, prepare kernels, and get kernel functions
+        self._set_macros()
+        self._set_field_knls()
+        if self.grid.pmls["slabs"]:
+            self._set_pml_knls()
+        if self.grid.rxs:
+            self._set_rx_knl()
+        if self.grid.voltagesources + self.grid.hertziandipoles + self.grid.magneticdipoles:
+            self._set_src_knls()
+        if self.grid.snapshots:
+            self._set_snapshot_knl()
+
+    def _build_knl(self, knl_func, subs_name_args, subs_func):
+        """Builds a CUDA kernel from templates: 1) function name and args;
+            and 2) function (kernel) body.
+
+        Args:
+            knl_func: dict containing templates for function name and args,
+                        and function body.
+            subs_name_args: dict containing substitutions to be used with
+                                function name and args.
+            subs_func: dict containing substitutions to be used with function
+                        (kernel) body.
+
+        Returns:
+            knl: string with complete kernel
+        """
+
+        name_plus_args = knl_func["args_cuda"].substitute(subs_name_args)
+        func_body = knl_func["func"].substitute(subs_func)
+        knl = self.knl_common + "\n" + name_plus_args + "{" + func_body + "}"
+
+        return knl
+
+    def _set_macros(self):
+        """Common macros to be used in kernels."""
+
+        # Set specific values for any dispersive materials
+        if config.get_model_config().materials["maxpoles"] > 0:
+            NY_MATDISPCOEFFS = self.grid.updatecoeffsdispersive.shape[1]
+            NX_T = self.grid.Tx.shape[1]
+            NY_T = self.grid.Tx.shape[2]
+            NZ_T = self.grid.Tx.shape[3]
+        else:  # Set to one any substitutions for dispersive materials.
+            NY_MATDISPCOEFFS = 1
+            NX_T = 1
+            NY_T = 1
+            NZ_T = 1
+
+        self.knl_common = self.env.get_template("knl_common_cuda.tmpl").render(
+            REAL=config.sim_config.dtypes["C_float_or_double"],
+            N_updatecoeffsE=self.grid.updatecoeffsE.size,
+            N_updatecoeffsH=self.grid.updatecoeffsH.size,
+            NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1],
+            NY_MATDISPCOEFFS=NY_MATDISPCOEFFS,
+            NX_FIELDS=self.grid.nx + 1,
+            NY_FIELDS=self.grid.ny + 1,
+            NZ_FIELDS=self.grid.nz + 1,
+            NX_ID=self.grid.ID.shape[1],
+            NY_ID=self.grid.ID.shape[2],
+            NZ_ID=self.grid.ID.shape[3],
+            NX_T=NX_T,
+            NY_T=NY_T,
+            NZ_T=NZ_T,
+            NY_RXCOORDS=3,
+            NX_RXS=6,
+            NY_RXS=self.grid.iterations,
+            NZ_RXS=len(self.grid.rxs),
+            NY_SRCINFO=4,
+            NY_SRCWAVES=self.grid.iterations,
+            NX_SNAPS=Snapshot.nx_max,
+            NY_SNAPS=Snapshot.ny_max,
+            NZ_SNAPS=Snapshot.nz_max,
+        )
+
+    def _set_field_knls(self):
+        """Electric and magnetic field updates - prepares kernels, and
+        gets kernel functions.
+        """
+
+        bld = self._build_knl(
+            knl_fields_updates.update_electric, self.subs_name_args, self.subs_func
+        )
+        knlE = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+        self.update_electric_dev = knlE.get_function("update_electric")
+
+        bld = self._build_knl(
+            knl_fields_updates.update_magnetic, self.subs_name_args, self.subs_func
+        )
+        knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+        self.update_magnetic_dev = knlH.get_function("update_magnetic")
+
+        self._copy_mat_coeffs(knlE, knlH)
+
+        # If there are any dispersive materials (updates are split into two
+        # parts as they require present and updated electric field values).
+        if config.get_model_config().materials["maxpoles"] > 0:
+            self.subs_func.update(
+                {
+                    "REAL": config.sim_config.dtypes["C_float_or_double"],
+                    "REALFUNC": config.get_model_config().materials["crealfunc"],
+                    "NX_T": self.grid.Tx.shape[1],
+                    "NY_T": self.grid.Tx.shape[2],
+                    "NZ_T": self.grid.Tx.shape[3],
+                }
+            )
+
+            bld = self._build_knl(
+                knl_fields_updates.update_electric_dispersive_A, self.subs_name_args, self.subs_func
+            )
+            knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            self.dispersive_update_a = knl.get_function("update_electric_dispersive_A")
+            self._copy_mat_coeffs(knl, knl)
+
+            bld = self._build_knl(
+                knl_fields_updates.update_electric_dispersive_B, self.subs_name_args, self.subs_func
+            )
+            knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            self.dispersive_update_b = knl.get_function("update_electric_dispersive_B")
+            self._copy_mat_coeffs(knl, knl)
+
+        # Set blocks per grid and initialise field arrays on GPU
+        self.grid.set_blocks_per_grid()
+        self.grid.htod_geometry_arrays()
+        self.grid.htod_field_arrays()
+        if config.get_model_config().materials["maxpoles"] > 0:
+            self.grid.htod_dispersive_arrays()
+
+    def _set_pml_knls(self):
+        """PMLS - prepares kernels and gets kernel functions."""
+        knl_pml_updates_electric = import_module(
+            "gprMax.cuda_opencl.knl_pml_updates_electric_" + self.grid.pmls["formulation"]
+        )
+        knl_pml_updates_magnetic = import_module(
+            "gprMax.cuda_opencl.knl_pml_updates_magnetic_" + self.grid.pmls["formulation"]
+        )
+
+        # Initialise arrays on GPU, set block per grid, and get kernel functions
+        for pml in self.grid.pmls["slabs"]:
+            pml.htod_field_arrays()
+            pml.set_blocks_per_grid()
+            knl_name = f"order{len(pml.CFS)}_{pml.direction}"
+            self.subs_name_args["FUNC"] = knl_name
+
+            knl_electric = getattr(knl_pml_updates_electric, knl_name)
+            bld = self._build_knl(knl_electric, self.subs_name_args, self.subs_func)
+            knlE = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            pml.update_electric_dev = knlE.get_function(knl_name)
+
+            knl_magnetic = getattr(knl_pml_updates_magnetic, knl_name)
+            bld = self._build_knl(knl_magnetic, self.subs_name_args, self.subs_func)
+            knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            pml.update_magnetic_dev = knlH.get_function(knl_name)
+
+            # Copy material coefficient arrays to constant memory of GPU - must
+            # be done for each kernel
+            self._copy_mat_coeffs(knlE, knlH)
+
+    def _set_rx_knl(self):
+        """Receivers - initialises arrays on GPU, prepares kernel and gets kernel
+        function.
+        """
+        self.rxcoords_dev, self.rxs_dev = htod_rx_arrays(self.grid)
+
+        self.subs_func.update(
+            {
+                "REAL": config.sim_config.dtypes["C_float_or_double"],
+                "NY_RXCOORDS": 3,
+                "NX_RXS": 6,
+                "NY_RXS": self.grid.iterations,
+                "NZ_RXS": len(self.grid.rxs),
+            }
+        )
+
+        bld = self._build_knl(knl_store_outputs.store_outputs, self.subs_name_args, self.subs_func)
+        knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+        self.store_outputs_dev = knl.get_function("store_outputs")
+
+    def _set_src_knls(self):
+        """Sources - initialises arrays on GPU, prepares kernel and gets kernel
+        function.
+        """
+        self.subs_func.update({"NY_SRCINFO": 4, "NY_SRCWAVES": self.grid.iteration})
+
+        if self.grid.hertziandipoles:
+            (
+                self.srcinfo1_hertzian_dev,
+                self.srcinfo2_hertzian_dev,
+                self.srcwaves_hertzian_dev,
+            ) = htod_src_arrays(self.grid.hertziandipoles, self.grid)
+            bld = self._build_knl(
+                knl_source_updates.update_hertzian_dipole, self.subs_name_args, self.subs_func
+            )
+            knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            self.update_hertzian_dipole_dev = knl.get_function("update_hertzian_dipole")
+        if self.grid.magneticdipoles:
+            (
+                self.srcinfo1_magnetic_dev,
+                self.srcinfo2_magnetic_dev,
+                self.srcwaves_magnetic_dev,
+            ) = htod_src_arrays(self.grid.magneticdipoles, self.grid)
+            bld = self._build_knl(
+                knl_source_updates.update_magnetic_dipole, self.subs_name_args, self.subs_func
+            )
+            knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            self.update_magnetic_dipole_dev = knl.get_function("update_magnetic_dipole")
+        if self.grid.voltagesources:
+            (
+                self.srcinfo1_voltage_dev,
+                self.srcinfo2_voltage_dev,
+                self.srcwaves_voltage_dev,
+            ) = htod_src_arrays(self.grid.voltagesources, self.grid)
+            bld = self._build_knl(
+                knl_source_updates.update_voltage_source, self.subs_name_args, self.subs_func
+            )
+            knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+            self.update_voltage_source_dev = knl.get_function("update_voltage_source")
+
+        self._copy_mat_coeffs(knl, knl)
+
+    def _set_snapshot_knl(self):
+        """Snapshots - initialises arrays on GPU, prepares kernel and gets kernel
+        function.
+        """
+        (
+            self.snapEx_dev,
+            self.snapEy_dev,
+            self.snapEz_dev,
+            self.snapHx_dev,
+            self.snapHy_dev,
+            self.snapHz_dev,
+        ) = htod_snapshot_array(self.grid)
+
+        self.subs_func.update(
+            {
+                "REAL": config.sim_config.dtypes["C_float_or_double"],
+                "NX_SNAPS": Snapshot.nx_max,
+                "NY_SNAPS": Snapshot.ny_max,
+                "NZ_SNAPS": Snapshot.nz_max,
+            }
+        )
+
+        bld = self._build_knl(knl_snapshots.store_snapshot, self.subs_name_args, self.subs_func)
+        knl = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"])
+        self.store_snapshot_dev = knl.get_function("store_snapshot")
+
+    def _copy_mat_coeffs(self, knlE, knlH):
+        """Copies material coefficient arrays to constant memory of GPU
+            (must be <64KB).
+
+        Args:
+            knlE: kernel for electric field.
+            knlH: kernel for magnetic field.
+        """
+
+        # Check if coefficient arrays will fit on constant memory of GPU
+        if (
+            self.grid.updatecoeffsE.nbytes + self.grid.updatecoeffsH.nbytes
+            > config.get_model_config().device["dev"].total_constant_memory
+        ):
+            device = config.get_model_config().device["dev"]
+            logger.exception(
+                f"Too many materials in the model to fit onto "
+                + f"constant memory of size {humanize.naturalsize(device.total_constant_memory)} "
+                + f"on {device.deviceID}: {' '.join(device.name().split())}"
+            )
+            raise ValueError
+
+        updatecoeffsE = knlE.get_global("updatecoeffsE")[0]
+        updatecoeffsH = knlH.get_global("updatecoeffsH")[0]
+        self.drv.memcpy_htod(updatecoeffsE, self.grid.updatecoeffsE)
+        self.drv.memcpy_htod(updatecoeffsH, self.grid.updatecoeffsH)
+
+    def store_outputs(self):
+        """Stores field component values for every receiver."""
+        if self.grid.rxs:
+            self.store_outputs_dev(
+                np.int32(len(self.grid.rxs)),
+                np.int32(self.grid.iteration),
+                self.rxcoords_dev.gpudata,
+                self.rxs_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                self.grid.Hx_dev.gpudata,
+                self.grid.Hy_dev.gpudata,
+                self.grid.Hz_dev.gpudata,
+                block=(1, 1, 1),
+                grid=(round32(len(self.grid.rxs)), 1, 1),
+            )
+
+    def store_snapshots(self, iteration):
+        """Stores any snapshots.
+
+        Args:
+            iteration: int for iteration number.
+        """
+
+        for i, snap in enumerate(self.grid.snapshots):
+            if snap.time == iteration + 1:
+                snapno = 0 if config.get_model_config().device["snapsgpu2cpu"] else i
+                self.store_snapshot_dev(
+                    np.int32(snapno),
+                    np.int32(snap.xs),
+                    np.int32(snap.xf),
+                    np.int32(snap.ys),
+                    np.int32(snap.yf),
+                    np.int32(snap.zs),
+                    np.int32(snap.zf),
+                    np.int32(snap.dx),
+                    np.int32(snap.dy),
+                    np.int32(snap.dz),
+                    self.grid.Ex_dev.gpudata,
+                    self.grid.Ey_dev.gpudata,
+                    self.grid.Ez_dev.gpudata,
+                    self.grid.Hx_dev.gpudata,
+                    self.grid.Hy_dev.gpudata,
+                    self.grid.Hz_dev.gpudata,
+                    self.snapEx_dev.gpudata,
+                    self.snapEy_dev.gpudata,
+                    self.snapEz_dev.gpudata,
+                    self.snapHx_dev.gpudata,
+                    self.snapHy_dev.gpudata,
+                    self.snapHz_dev.gpudata,
+                    block=Snapshot.tpb,
+                    grid=Snapshot.bpg,
+                )
+                if config.get_model_config().device["snapsgpu2cpu"]:
+                    dtoh_snapshot_array(
+                        self.snapEx_dev.get(),
+                        self.snapEy_dev.get(),
+                        self.snapEz_dev.get(),
+                        self.snapHx_dev.get(),
+                        self.snapHy_dev.get(),
+                        self.snapHz_dev.get(),
+                        0,
+                        snap,
+                    )
+
+    def update_magnetic(self):
+        """Updates magnetic field components."""
+        self.update_magnetic_dev(
+            np.int32(self.grid.nx),
+            np.int32(self.grid.ny),
+            np.int32(self.grid.nz),
+            self.grid.ID_dev.gpudata,
+            self.grid.Hx_dev.gpudata,
+            self.grid.Hy_dev.gpudata,
+            self.grid.Hz_dev.gpudata,
+            self.grid.Ex_dev.gpudata,
+            self.grid.Ey_dev.gpudata,
+            self.grid.Ez_dev.gpudata,
+            block=self.grid.tpb,
+            grid=self.grid.bpg,
+        )
+
+    def update_magnetic_pml(self):
+        """Updates magnetic field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_magnetic()
+
+    def update_magnetic_sources(self):
+        """Updates magnetic field components from sources."""
+        if self.grid.magneticdipoles:
+            self.update_magnetic_dipole_dev(
+                np.int32(len(self.grid.magneticdipoles)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_magnetic_dev.gpudata,
+                self.srcinfo2_magnetic_dev.gpudata,
+                self.srcwaves_magnetic_dev.gpudata,
+                self.grid.ID_dev.gpudata,
+                self.grid.Hx_dev.gpudata,
+                self.grid.Hy_dev.gpudata,
+                self.grid.Hz_dev.gpudata,
+                block=(1, 1, 1),
+                grid=(round32(len(self.grid.magneticdipoles)), 1, 1),
+            )
+
+    def update_electric_a(self):
+        """Updates electric field components."""
+        # All materials are non-dispersive so do standard update.
+        if config.get_model_config().materials["maxpoles"] == 0:
+            self.update_electric_dev(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                self.grid.ID_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                self.grid.Hx_dev.gpudata,
+                self.grid.Hy_dev.gpudata,
+                self.grid.Hz_dev.gpudata,
+                block=self.grid.tpb,
+                grid=self.grid.bpg,
+            )
+
+        # 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).
+        else:
+            self.dispersive_update_a(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                np.int32(config.get_model_config().materials["maxpoles"]),
+                self.grid.updatecoeffsdispersive_dev.gpudata,
+                self.grid.Tx_dev.gpudata,
+                self.grid.Ty_dev.gpudata,
+                self.grid.Tz_dev.gpudata,
+                self.grid.ID_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                self.grid.Hx_dev.gpudata,
+                self.grid.Hy_dev.gpudata,
+                self.grid.Hz_dev.gpudata,
+                block=self.grid.tpb,
+                grid=self.grid.bpg,
+            )
+
+    def update_electric_pml(self):
+        """Updates electric field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_electric()
+
+    def update_electric_sources(self):
+        """Updates electric field components from sources -
+        update any Hertzian dipole sources last.
+        """
+        if self.grid.voltagesources:
+            self.update_voltage_source_dev(
+                np.int32(len(self.grid.voltagesources)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_voltage_dev.gpudata,
+                self.srcinfo2_voltage_dev.gpudata,
+                self.srcwaves_voltage_dev.gpudata,
+                self.grid.ID_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                block=(1, 1, 1),
+                grid=(round32(len(self.grid.voltagesources)), 1, 1),
+            )
+
+        if self.grid.hertziandipoles:
+            self.update_hertzian_dipole_dev(
+                np.int32(len(self.grid.hertziandipoles)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_hertzian_dev.gpudata,
+                self.srcinfo2_hertzian_dev.gpudata,
+                self.srcwaves_hertzian_dev.gpudata,
+                self.grid.ID_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                block=(1, 1, 1),
+                grid=(round32(len(self.grid.hertziandipoles)), 1, 1),
+            )
+
+        self.grid.iteration += 1
+
+    def update_electric_b(self):
+        """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 config.get_model_config().materials["maxpoles"] > 0:
+            self.dispersive_update_b(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                np.int32(config.get_model_config().materials["maxpoles"]),
+                self.grid.updatecoeffsdispersive_dev.gpudata,
+                self.grid.Tx_dev.gpudata,
+                self.grid.Ty_dev.gpudata,
+                self.grid.Tz_dev.gpudata,
+                self.grid.ID_dev.gpudata,
+                self.grid.Ex_dev.gpudata,
+                self.grid.Ey_dev.gpudata,
+                self.grid.Ez_dev.gpudata,
+                block=self.grid.tpb,
+                grid=self.grid.bpg,
+            )
+
+    def time_start(self):
+        """Starts event timers used to calculate solving time for model."""
+        self.iterstart = self.drv.Event()
+        self.iterend = self.drv.Event()
+        self.iterstart.record()
+        self.iterstart.synchronize()
+
+    def calculate_memory_used(self, iteration):
+        """Calculates memory used on last iteration.
+
+        Args:
+            iteration: int for iteration number.
+
+        Returns:
+            Memory (RAM) used on GPU.
+        """
+        if iteration == self.grid.iterations - 1:
+            # Total minus free memory in current context
+            return self.drv.mem_get_info()[1] - self.drv.mem_get_info()[0]
+
+    def calculate_solve_time(self):
+        """Calculates solving time for model."""
+        self.iterend.record()
+        self.iterend.synchronize()
+        return self.iterstart.time_till(self.iterend) * 1e-3
+
+    def finalise(self):
+        """Copies data from GPU back to CPU to save to file(s)."""
+        # Copy output from receivers array back to correct receiver objects
+        if self.grid.rxs:
+            dtoh_rx_array(self.rxs_dev.get(), self.rxcoords_dev.get(), self.grid)
+
+        # Copy data from any snapshots back to correct snapshot objects
+        if self.grid.snapshots and not config.get_model_config().device["snapsgpu2cpu"]:
+            for i, snap in enumerate(self.grid.snapshots):
+                dtoh_snapshot_array(
+                    self.snapEx_dev.get(),
+                    self.snapEy_dev.get(),
+                    self.snapEz_dev.get(),
+                    self.snapHx_dev.get(),
+                    self.snapHy_dev.get(),
+                    self.snapHz_dev.get(),
+                    i,
+                    snap,
+                )
+
+    def cleanup(self):
+        """Cleanup GPU context."""
+        # Remove context from top of stack and clear
+        self.ctx.pop()
+        self.ctx = None

gprMax/updates/opencl_updates.py

@@ -0,0 +1,615 @@
+# Copyright (C) 2015-2024: 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 importlib import import_module
+
+import numpy as np
+from jinja2 import Environment, PackageLoader
+
+from gprMax import config
+from gprMax.cuda_opencl import (
+    knl_fields_updates,
+    knl_snapshots,
+    knl_source_updates,
+    knl_store_outputs,
+)
+from gprMax.grid.opencl_grid import OpenCLGrid
+from gprMax.receivers import dtoh_rx_array, htod_rx_arrays
+from gprMax.snapshots import Snapshot, dtoh_snapshot_array, htod_snapshot_array
+from gprMax.sources import htod_src_arrays
+from gprMax.updates.updates import Updates
+
+logger = logging.getLogger(__name__)
+
+
+class OpenCLUpdates(Updates):
+    """Defines update functions for OpenCL-based solver."""
+
+    def __init__(self, G: OpenCLGrid):
+        """
+        Args:
+            G: OpenCLGrid class describing a grid in a model.
+        """
+
+        self.grid = G
+
+        # Import pyopencl module
+        self.cl = import_module("pyopencl")
+        self.elwiseknl = getattr(import_module("pyopencl.elementwise"), "ElementwiseKernel")
+
+        # Select device, create context and command queue
+        self.dev = config.get_model_config().device["dev"]
+        self.ctx = self.cl.Context(devices=[self.dev])
+        self.queue = self.cl.CommandQueue(
+            self.ctx, properties=self.cl.command_queue_properties.PROFILING_ENABLE
+        )
+
+        # Enviroment for templating kernels
+        self.env = Environment(loader=PackageLoader("gprMax", "cuda_opencl"))
+
+        # Initialise arrays on device, prepare kernels, and get kernel functions
+        self._set_macros()
+        self._set_field_knls()
+        if self.grid.pmls["slabs"]:
+            self._set_pml_knls()
+        if self.grid.rxs:
+            self._set_rx_knl()
+        if self.grid.voltagesources + self.grid.hertziandipoles + self.grid.magneticdipoles:
+            self._set_src_knls()
+        if self.grid.snapshots:
+            self._set_snapshot_knl()
+
+    def _set_macros(self):
+        """Common macros to be used in kernels."""
+
+        # Set specific values for any dispersive materials
+        if config.get_model_config().materials["maxpoles"] > 0:
+            NY_MATDISPCOEFFS = self.grid.updatecoeffsdispersive.shape[1]
+            NX_T = self.grid.Tx.shape[1]
+            NY_T = self.grid.Tx.shape[2]
+            NZ_T = self.grid.Tx.shape[3]
+        else:  # Set to one any substitutions for dispersive materials.
+            NY_MATDISPCOEFFS = 1
+            NX_T = 1
+            NY_T = 1
+            NZ_T = 1
+
+        self.knl_common = self.env.get_template("knl_common_opencl.tmpl").render(
+            updatecoeffsE=self.grid.updatecoeffsE.ravel(),
+            updatecoeffsH=self.grid.updatecoeffsH.ravel(),
+            REAL=config.sim_config.dtypes["C_float_or_double"],
+            N_updatecoeffsE=self.grid.updatecoeffsE.size,
+            N_updatecoeffsH=self.grid.updatecoeffsH.size,
+            NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1],
+            NY_MATDISPCOEFFS=NY_MATDISPCOEFFS,
+            NX_FIELDS=self.grid.nx + 1,
+            NY_FIELDS=self.grid.ny + 1,
+            NZ_FIELDS=self.grid.nz + 1,
+            NX_ID=self.grid.ID.shape[1],
+            NY_ID=self.grid.ID.shape[2],
+            NZ_ID=self.grid.ID.shape[3],
+            NX_T=NX_T,
+            NY_T=NY_T,
+            NZ_T=NZ_T,
+            NY_RXCOORDS=3,
+            NX_RXS=6,
+            NY_RXS=self.grid.iterations,
+            NZ_RXS=len(self.grid.rxs),
+            NY_SRCINFO=4,
+            NY_SRCWAVES=self.grid.iterations,
+            NX_SNAPS=Snapshot.nx_max,
+            NY_SNAPS=Snapshot.ny_max,
+            NZ_SNAPS=Snapshot.nz_max,
+        )
+
+    def _set_field_knls(self):
+        """Electric and magnetic field updates - prepares kernels, and
+        gets kernel functions.
+        """
+
+        subs = {
+            "CUDA_IDX": "",
+            "NX_FIELDS": self.grid.nx + 1,
+            "NY_FIELDS": self.grid.ny + 1,
+            "NZ_FIELDS": self.grid.nz + 1,
+            "NX_ID": self.grid.ID.shape[1],
+            "NY_ID": self.grid.ID.shape[2],
+            "NZ_ID": self.grid.ID.shape[3],
+        }
+
+        self.update_electric_dev = self.elwiseknl(
+            self.ctx,
+            knl_fields_updates.update_electric["args_opencl"].substitute(
+                {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+            ),
+            knl_fields_updates.update_electric["func"].substitute(subs),
+            "update_electric",
+            preamble=self.knl_common,
+            options=config.sim_config.devices["compiler_opts"],
+        )
+
+        self.update_magnetic_dev = self.elwiseknl(
+            self.ctx,
+            knl_fields_updates.update_magnetic["args_opencl"].substitute(
+                {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+            ),
+            knl_fields_updates.update_magnetic["func"].substitute(subs),
+            "update_magnetic",
+            preamble=self.knl_common,
+            options=config.sim_config.devices["compiler_opts"],
+        )
+
+        # If there are any dispersive materials (updates are split into two
+        # parts as they require present and updated electric field values).
+        if config.get_model_config().materials["maxpoles"] > 0:
+            subs = {
+                "CUDA_IDX": "",
+                "REAL": config.sim_config.dtypes["C_float_or_double"],
+                "REALFUNC": config.get_model_config().materials["crealfunc"],
+                "NX_FIELDS": self.grid.nx + 1,
+                "NY_FIELDS": self.grid.ny + 1,
+                "NZ_FIELDS": self.grid.nz + 1,
+                "NX_ID": self.grid.ID.shape[1],
+                "NY_ID": self.grid.ID.shape[2],
+                "NZ_ID": self.grid.ID.shape[3],
+                "NX_T": self.grid.Tx.shape[1],
+                "NY_T": self.grid.Tx.shape[2],
+                "NZ_T": self.grid.Tx.shape[3],
+            }
+
+            self.dispersive_update_a = self.elwiseknl(
+                self.ctx,
+                knl_fields_updates.update_electric_dispersive_A["args_opencl"].substitute(
+                    {
+                        "REAL": config.sim_config.dtypes["C_float_or_double"],
+                        "COMPLEX": config.get_model_config().materials["dispersiveCdtype"],
+                    }
+                ),
+                knl_fields_updates.update_electric_dispersive_A["func"].substitute(subs),
+                "update_electric_dispersive_A",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+
+            self.dispersive_update_b = self.elwiseknl(
+                self.ctx,
+                knl_fields_updates.update_electric_dispersive_B["args_opencl"].substitute(
+                    {
+                        "REAL": config.sim_config.dtypes["C_float_or_double"],
+                        "COMPLEX": config.get_model_config().materials["dispersiveCdtype"],
+                    }
+                ),
+                knl_fields_updates.update_electric_dispersive_B["func"].substitute(subs),
+                "update_electric_dispersive_B",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+
+        # Initialise field arrays on compute device
+        self.grid.htod_geometry_arrays(self.queue)
+        self.grid.htod_field_arrays(self.queue)
+        if config.get_model_config().materials["maxpoles"] > 0:
+            self.grid.htod_dispersive_arrays(self.queue)
+
+    def _set_pml_knls(self):
+        """PMLS - prepares kernels and gets kernel functions."""
+        knl_pml_updates_electric = import_module(
+            "gprMax.cuda_opencl.knl_pml_updates_electric_" + self.grid.pmls["formulation"]
+        )
+        knl_pml_updates_magnetic = import_module(
+            "gprMax.cuda_opencl.knl_pml_updates_magnetic_" + self.grid.pmls["formulation"]
+        )
+
+        subs = {
+            "CUDA_IDX": "",
+            "REAL": config.sim_config.dtypes["C_float_or_double"],
+            "NX_FIELDS": self.grid.nx + 1,
+            "NY_FIELDS": self.grid.ny + 1,
+            "NZ_FIELDS": self.grid.nz + 1,
+            "NX_ID": self.grid.ID.shape[1],
+            "NY_ID": self.grid.ID.shape[2],
+            "NZ_ID": self.grid.ID.shape[3],
+        }
+
+        # Set workgroup size, initialise arrays on compute device, and get
+        # kernel functions
+        for pml in self.grid.pmls["slabs"]:
+            pml.set_queue(self.queue)
+            pml.htod_field_arrays()
+            knl_name = f"order{len(pml.CFS)}_{pml.direction}"
+            knl_electric_name = getattr(knl_pml_updates_electric, knl_name)
+            knl_magnetic_name = getattr(knl_pml_updates_magnetic, knl_name)
+
+            pml.update_electric_dev = self.elwiseknl(
+                self.ctx,
+                knl_electric_name["args_opencl"].substitute(
+                    {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                knl_electric_name["func"].substitute(subs),
+                f"pml_updates_electric_{knl_name}",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+
+            pml.update_magnetic_dev = self.elwiseknl(
+                self.ctx,
+                knl_magnetic_name["args_opencl"].substitute(
+                    {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                knl_magnetic_name["func"].substitute(subs),
+                f"pml_updates_magnetic_{knl_name}",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+
+    def _set_rx_knl(self):
+        """Receivers - initialises arrays on compute device, prepares kernel and
+        gets kernel function.
+        """
+        self.rxcoords_dev, self.rxs_dev = htod_rx_arrays(self.grid, self.queue)
+        self.store_outputs_dev = self.elwiseknl(
+            self.ctx,
+            knl_store_outputs.store_outputs["args_opencl"].substitute(
+                {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+            ),
+            knl_store_outputs.store_outputs["func"].substitute({"CUDA_IDX": ""}),
+            "store_outputs",
+            preamble=self.knl_common,
+            options=config.sim_config.devices["compiler_opts"],
+        )
+
+    def _set_src_knls(self):
+        """Sources - initialises arrays on compute device, prepares kernel and
+        gets kernel function.
+        """
+        if self.grid.hertziandipoles:
+            (
+                self.srcinfo1_hertzian_dev,
+                self.srcinfo2_hertzian_dev,
+                self.srcwaves_hertzian_dev,
+            ) = htod_src_arrays(self.grid.hertziandipoles, self.grid, self.queue)
+            self.update_hertzian_dipole_dev = self.elwiseknl(
+                self.ctx,
+                knl_source_updates.update_hertzian_dipole["args_opencl"].substitute(
+                    {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                knl_source_updates.update_hertzian_dipole["func"].substitute(
+                    {"CUDA_IDX": "", "REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                "update_hertzian_dipole",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+        if self.grid.magneticdipoles:
+            (
+                self.srcinfo1_magnetic_dev,
+                self.srcinfo2_magnetic_dev,
+                self.srcwaves_magnetic_dev,
+            ) = htod_src_arrays(self.grid.magneticdipoles, self.grid, self.queue)
+            self.update_magnetic_dipole_dev = self.elwiseknl(
+                self.ctx,
+                knl_source_updates.update_magnetic_dipole["args_opencl"].substitute(
+                    {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                knl_source_updates.update_magnetic_dipole["func"].substitute(
+                    {"CUDA_IDX": "", "REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                "update_magnetic_dipole",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+        if self.grid.voltagesources:
+            (
+                self.srcinfo1_voltage_dev,
+                self.srcinfo2_voltage_dev,
+                self.srcwaves_voltage_dev,
+            ) = htod_src_arrays(self.grid.voltagesources, self.grid, self.queue)
+            self.update_voltage_source_dev = self.elwiseknl(
+                self.ctx,
+                knl_source_updates.update_voltage_source["args_opencl"].substitute(
+                    {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                knl_source_updates.update_voltage_source["func"].substitute(
+                    {"CUDA_IDX": "", "REAL": config.sim_config.dtypes["C_float_or_double"]}
+                ),
+                "update_voltage_source",
+                preamble=self.knl_common,
+                options=config.sim_config.devices["compiler_opts"],
+            )
+
+    def _set_snapshot_knl(self):
+        """Snapshots - initialises arrays on compute device, prepares kernel and
+        gets kernel function.
+        """
+        (
+            self.snapEx_dev,
+            self.snapEy_dev,
+            self.snapEz_dev,
+            self.snapHx_dev,
+            self.snapHy_dev,
+            self.snapHz_dev,
+        ) = htod_snapshot_array(self.grid, self.queue)
+        self.store_snapshot_dev = self.elwiseknl(
+            self.ctx,
+            knl_snapshots.store_snapshot["args_opencl"].substitute(
+                {"REAL": config.sim_config.dtypes["C_float_or_double"]}
+            ),
+            knl_snapshots.store_snapshot["func"].substitute(
+                {
+                    "CUDA_IDX": "",
+                    "NX_SNAPS": Snapshot.nx_max,
+                    "NY_SNAPS": Snapshot.ny_max,
+                    "NZ_SNAPS": Snapshot.nz_max,
+                }
+            ),
+            "store_snapshot",
+            preamble=self.knl_common,
+            options=config.sim_config.devices["compiler_opts"],
+        )
+
+    def store_outputs(self):
+        """Stores field component values for every receiver."""
+        if self.grid.rxs:
+            self.store_outputs_dev(
+                np.int32(len(self.grid.rxs)),
+                np.int32(self.grid.iteration),
+                self.rxcoords_dev,
+                self.rxs_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+                self.grid.Hx_dev,
+                self.grid.Hy_dev,
+                self.grid.Hz_dev,
+            )
+
+    def store_snapshots(self, iteration):
+        """Stores any snapshots.
+
+        Args:
+            iteration: int for iteration number.
+        """
+
+        for i, snap in enumerate(self.grid.snapshots):
+            if snap.time == iteration + 1:
+                snapno = 0 if config.get_model_config().device["snapsgpu2cpu"] else i
+                self.store_snapshot_dev(
+                    np.int32(snapno),
+                    np.int32(snap.xs),
+                    np.int32(snap.xf),
+                    np.int32(snap.ys),
+                    np.int32(snap.yf),
+                    np.int32(snap.zs),
+                    np.int32(snap.zf),
+                    np.int32(snap.dx),
+                    np.int32(snap.dy),
+                    np.int32(snap.dz),
+                    self.grid.Ex_dev,
+                    self.grid.Ey_dev,
+                    self.grid.Ez_dev,
+                    self.grid.Hx_dev,
+                    self.grid.Hy_dev,
+                    self.grid.Hz_dev,
+                    self.snapEx_dev,
+                    self.snapEy_dev,
+                    self.snapEz_dev,
+                    self.snapHx_dev,
+                    self.snapHy_dev,
+                    self.snapHz_dev,
+                )
+
+                if config.get_model_config().device["snapsgpu2cpu"]:
+                    dtoh_snapshot_array(
+                        self.snapEx_dev.get(),
+                        self.snapEy_dev.get(),
+                        self.snapEz_dev.get(),
+                        self.snapHx_dev.get(),
+                        self.snapHy_dev.get(),
+                        self.snapHz_dev.get(),
+                        0,
+                        snap,
+                    )
+
+    def update_magnetic(self):
+        """Updates magnetic field components."""
+        self.update_magnetic_dev(
+            np.int32(self.grid.nx),
+            np.int32(self.grid.ny),
+            np.int32(self.grid.nz),
+            self.grid.ID_dev,
+            self.grid.Hx_dev,
+            self.grid.Hy_dev,
+            self.grid.Hz_dev,
+            self.grid.Ex_dev,
+            self.grid.Ey_dev,
+            self.grid.Ez_dev,
+        )
+
+    def update_magnetic_pml(self):
+        """Updates magnetic field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_magnetic()
+
+    def update_magnetic_sources(self):
+        """Updates magnetic field components from sources."""
+        if self.grid.magneticdipoles:
+            self.update_magnetic_dipole_dev(
+                np.int32(len(self.grid.magneticdipoles)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_magnetic_dev,
+                self.srcinfo2_magnetic_dev,
+                self.srcwaves_magnetic_dev,
+                self.grid.ID_dev,
+                self.grid.Hx_dev,
+                self.grid.Hy_dev,
+                self.grid.Hz_dev,
+            )
+
+    def update_electric_a(self):
+        """Updates electric field components."""
+        # All materials are non-dispersive so do standard update.
+        if config.get_model_config().materials["maxpoles"] == 0:
+            self.update_electric_dev(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                self.grid.ID_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+                self.grid.Hx_dev,
+                self.grid.Hy_dev,
+                self.grid.Hz_dev,
+            )
+
+        # 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).
+        else:
+            self.dispersive_update_a(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                np.int32(config.get_model_config().materials["maxpoles"]),
+                self.grid.ID_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+                self.grid.Hx_dev,
+                self.grid.Hy_dev,
+                self.grid.Hz_dev,
+                self.grid.updatecoeffsdispersive_dev,
+                self.grid.Tx_dev,
+                self.grid.Ty_dev,
+                self.grid.Tz_dev,
+            )
+
+    def update_electric_pml(self):
+        """Updates electric field components with the PML correction."""
+        for pml in self.grid.pmls["slabs"]:
+            pml.update_electric()
+
+    def update_electric_sources(self):
+        """Updates electric field components from sources -
+        update any Hertzian dipole sources last.
+        """
+        if self.grid.voltagesources:
+            self.update_voltage_source_dev(
+                np.int32(len(self.grid.voltagesources)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_voltage_dev,
+                self.srcinfo2_voltage_dev,
+                self.srcwaves_voltage_dev,
+                self.grid.ID_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+            )
+
+        if self.grid.hertziandipoles:
+            self.update_hertzian_dipole_dev(
+                np.int32(len(self.grid.hertziandipoles)),
+                np.int32(self.grid.iteration),
+                config.sim_config.dtypes["float_or_double"](self.grid.dx),
+                config.sim_config.dtypes["float_or_double"](self.grid.dy),
+                config.sim_config.dtypes["float_or_double"](self.grid.dz),
+                self.srcinfo1_hertzian_dev,
+                self.srcinfo2_hertzian_dev,
+                self.srcwaves_hertzian_dev,
+                self.grid.ID_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+            )
+
+        self.grid.iteration += 1
+
+    def update_electric_b(self):
+        """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 config.get_model_config().materials["maxpoles"] > 0:
+            self.dispersive_update_b(
+                np.int32(self.grid.nx),
+                np.int32(self.grid.ny),
+                np.int32(self.grid.nz),
+                np.int32(config.get_model_config().materials["maxpoles"]),
+                self.grid.ID_dev,
+                self.grid.Ex_dev,
+                self.grid.Ey_dev,
+                self.grid.Ez_dev,
+                self.grid.updatecoeffsdispersive_dev,
+                self.grid.Tx_dev,
+                self.grid.Ty_dev,
+                self.grid.Tz_dev,
+            )
+
+    def time_start(self):
+        """Starts event timers used to calculate solving time for model."""
+        self.event_marker1 = self.cl.enqueue_marker(self.queue)
+        self.event_marker1.wait()
+
+    def calculate_memory_used(self, iteration):
+        """Calculates memory used on last iteration.
+
+        Args:
+            iteration: int for iteration number.
+
+        Returns:
+            Memory (RAM) used on compute device.
+        """
+        # No clear way to determine memory used from PyOpenCL unlike PyCUDA.
+        pass
+
+    def calculate_solve_time(self):
+        """Calculates solving time for model."""
+        event_marker2 = self.cl.enqueue_marker(self.queue)
+        event_marker2.wait()
+        return (event_marker2.profile.end - self.event_marker1.profile.start) * 1e-9
+
+    def finalise(self):
+        """Copies data from compute device back to CPU to save to file(s)."""
+        # Copy output from receivers array back to correct receiver objects
+        if self.grid.rxs:
+            dtoh_rx_array(self.rxs_dev.get(), self.rxcoords_dev.get(), self.grid)
+
+        # Copy data from any snapshots back to correct snapshot objects
+        if self.grid.snapshots and not config.get_model_config().device["snapsgpu2cpu"]:
+            for i, snap in enumerate(self.grid.snapshots):
+                dtoh_snapshot_array(
+                    self.snapEx_dev.get(),
+                    self.snapEy_dev.get(),
+                    self.snapEz_dev.get(),
+                    self.snapHx_dev.get(),
+                    self.snapHy_dev.get(),
+                    self.snapHz_dev.get(),
+                    i,
+                    snap,
+                )
+
+    def cleanup(self):
+        pass

reframe_tests/job_scripts/archer2_tests.slurm

@@ -11,6 +11,8 @@
 #   using threading.
 export OMP_NUM_THREADS=1
 
+module load cray-python
+
 source ../.venv/bin/activate
 
 reframe -C configuration/archer2_settings.py -c reframe_tests.py -c base_tests.py -r "$@"

requirements.txt

@@ -12,7 +12,6 @@ matplotlib
 numpy
 numpy-stl
 pandas
-pip
 pre-commit
 psutil
 # pycuda

setup.py

@@ -253,7 +253,7 @@ else:
     )
 
     # Parse long_description from README.rst file.
-    with open("README.rst", "r") as fd:
+    with open("README.rst", "r", encoding="utf-8") as fd:
         long_description = fd.read()
 
     setup(

tests/test_update_functions.py

@@ -1,116 +0,0 @@
-import argparse
-
-import numpy as np
-import pytest
-
-from gprMax import config, gprMax
-from gprMax.grid import FDTDGrid
-from gprMax.materials import create_built_in_materials
-from gprMax.model_build_run import GridBuilder
-from gprMax.pml import CFS
-from gprMax.updates.updates import CPUUpdates
-
-
-def build_grid(nx, ny, nz, dl=0.001, dt=3e-9):
-    grid = FDTDGrid()
-    grid.nx = nx
-    grid.ny = ny
-    grid.nz = nz
-    grid.dx = dl
-    grid.dy = dl
-    grid.dz = dl
-    grid.dt = dt
-    create_built_in_materials(grid)
-    grid.initialise_geometry_arrays()
-    grid.initialise_field_arrays()
-    grid.initialise_std_update_coeff_arrays()
-    grid.pmls["cfs"] = [CFS()]
-
-    grid_builder = GridBuilder(grid)
-    grid_builder.build_pmls()
-    grid_builder.build_components()
-    grid_builder.build_materials()
-
-    return grid
-
-
-@pytest.fixture
-def config_mock(monkeypatch):
-    def _mock_simulation_config():
-        args = argparse.Namespace(**gprMax.args_defaults)
-        args.inputfile = "test.in"
-        return config.SimulationConfig(args)
-
-    def _mock_model_config():
-        model_config = config.ModelConfig()
-        model_config.ompthreads = 1
-        return model_config
-
-    monkeypatch.setattr(config, "sim_config", _mock_simulation_config())
-    monkeypatch.setattr(config, "get_model_config", _mock_model_config)
-
-
-def test_update_magnetic_cpu(config_mock):
-    grid = build_grid(100, 100, 100)
-
-    expected_value = np.zeros((101, 101, 101))
-
-    cpu_updates = CPUUpdates(grid)
-    cpu_updates.update_magnetic()
-
-    assert np.equal(grid.Ex, expected_value).all()
-    assert np.equal(grid.Ey, expected_value).all()
-    assert np.equal(grid.Ez, expected_value).all()
-    assert np.equal(grid.Hx, expected_value).all()
-    assert np.equal(grid.Hy, expected_value).all()
-    assert np.equal(grid.Hz, expected_value).all()
-
-    for pml in grid.pmls["slabs"]:
-        assert np.equal(pml.HPhi1, 0).all()
-        assert np.equal(pml.HPhi2, 0).all()
-        assert np.equal(pml.EPhi1, 0).all()
-        assert np.equal(pml.EPhi2, 0).all()
-
-
-def test_update_magnetic_pml_cpu(config_mock):
-    grid = build_grid(100, 100, 100)
-
-    grid_expected_value = np.zeros((101, 101, 101))
-
-    cpu_updates = CPUUpdates(grid)
-    cpu_updates.update_magnetic_pml()
-
-    assert np.equal(grid.Ex, grid_expected_value).all()
-    assert np.equal(grid.Ey, grid_expected_value).all()
-    assert np.equal(grid.Ez, grid_expected_value).all()
-    assert np.equal(grid.Hx, grid_expected_value).all()
-    assert np.equal(grid.Hy, grid_expected_value).all()
-    assert np.equal(grid.Hz, grid_expected_value).all()
-
-    for pml in grid.pmls["slabs"]:
-        assert np.equal(pml.HPhi1, 0).all()
-        assert np.equal(pml.HPhi2, 0).all()
-        assert np.equal(pml.EPhi1, 0).all()
-        assert np.equal(pml.EPhi2, 0).all()
-
-
-def test_update_electric_pml_cpu(config_mock):
-    grid = build_grid(100, 100, 100)
-
-    grid_expected_value = np.zeros((101, 101, 101))
-
-    cpu_updates = CPUUpdates(grid)
-    cpu_updates.update_electric_pml()
-
-    assert np.equal(grid.Ex, grid_expected_value).all()
-    assert np.equal(grid.Ey, grid_expected_value).all()
-    assert np.equal(grid.Ez, grid_expected_value).all()
-    assert np.equal(grid.Hx, grid_expected_value).all()
-    assert np.equal(grid.Hy, grid_expected_value).all()
-    assert np.equal(grid.Hz, grid_expected_value).all()
-
-    for pml in grid.pmls["slabs"]:
-        assert np.equal(pml.HPhi1, 0).all()
-        assert np.equal(pml.HPhi2, 0).all()
-        assert np.equal(pml.EPhi1, 0).all()
-        assert np.equal(pml.EPhi2, 0).all()

tests/updates/test_cpu_updates.py

@@ -0,0 +1,244 @@
+import argparse
+
+import numpy as np
+import pytest
+from pytest import MonkeyPatch
+
+from gprMax import config, gprMax
+from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.materials import create_built_in_materials
+from gprMax.model import GridBuilder
+from gprMax.pml import CFS
+from gprMax.updates.cpu_updates import CPUUpdates
+
+
+def build_grid(nx: int, ny: int, nz: int, dl: float = 0.001, dt: float = 3e-9) -> FDTDGrid:
+    grid = FDTDGrid()
+    grid.nx = nx
+    grid.ny = ny
+    grid.nz = nz
+    grid.dx = dl
+    grid.dy = dl
+    grid.dz = dl
+    grid.dt = dt
+    create_built_in_materials(grid)
+    grid.initialise_geometry_arrays()
+    grid.initialise_field_arrays()
+    grid.initialise_std_update_coeff_arrays()
+    grid.pmls["cfs"] = [CFS()]
+
+    grid_builder = GridBuilder(grid)
+    grid_builder.build_pmls()
+    grid_builder.build_components()
+    grid_builder.build_materials()
+
+    return grid
+
+
+@pytest.fixture
+def config_mock(monkeypatch: MonkeyPatch):
+    def _mock_simulation_config() -> config.SimulationConfig:
+        args = argparse.Namespace(**gprMax.args_defaults)
+        args.inputfile = "test.in"
+        return config.SimulationConfig(args)
+
+    def _mock_model_config() -> config.ModelConfig:
+        model_config = config.ModelConfig(1)
+        model_config.ompthreads = 1
+        return model_config
+
+    monkeypatch.setattr(config, "sim_config", _mock_simulation_config())
+    monkeypatch.setattr(config, "get_model_config", _mock_model_config)
+
+
+def test_update_magnetic_zero_grid(config_mock):
+    grid = build_grid(100, 100, 100)
+
+    expected_value = np.zeros((101, 101, 101))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_magnetic()
+
+    assert np.equal(grid.Ex, expected_value).all()
+    assert np.equal(grid.Ey, expected_value).all()
+    assert np.equal(grid.Ez, expected_value).all()
+    assert np.equal(grid.Hx, expected_value).all()
+    assert np.equal(grid.Hy, expected_value).all()
+    assert np.equal(grid.Hz, expected_value).all()
+
+    for pml in grid.pmls["slabs"]:
+        assert np.equal(pml.HPhi1, 0).all()
+        assert np.equal(pml.HPhi2, 0).all()
+        assert np.equal(pml.EPhi1, 0).all()
+        assert np.equal(pml.EPhi2, 0).all()
+
+
+def test_update_magnetic(config_mock):
+    grid = build_grid(11, 11, 11)
+
+    grid.updatecoeffsH[1] = 1
+
+    grid.Ex = np.tile(np.array([[[1, 2], [2, 1]], [[2, 1], [1, 2]]], dtype=np.float32), (6, 6, 6))
+    grid.Ey = np.tile(np.array([[[1, 3], [3, 1]], [[3, 1], [1, 3]]], dtype=np.float32), (6, 6, 6))
+    grid.Ez = np.tile(np.array([[[1, 4], [4, 1]], [[4, 1], [1, 4]]], dtype=np.float32), (6, 6, 6))
+    grid.Hx = np.tile(np.array([[[3, 1], [1, 3]], [[1, 3], [3, 1]]], dtype=np.float32), (6, 6, 6))
+    grid.Hy = np.tile(np.array([[[1, 5], [5, 1]], [[5, 1], [1, 5]]], dtype=np.float32), (6, 6, 6))
+    grid.Hz = np.tile(np.array([[[5, 3], [3, 5]], [[3, 5], [5, 3]]], dtype=np.float32), (6, 6, 6))
+
+    expected_Ex = grid.Ex.copy()
+    expected_Ey = grid.Ey.copy()
+    expected_Ez = grid.Ez.copy()
+    expected_Hx = grid.Hx.copy()
+    expected_Hy = grid.Hy.copy()
+    expected_Hz = grid.Hz.copy()
+    expected_Hx[1:, :-1, :-1] = np.tile(np.array([[[2]]], dtype=np.float32), (11, 11, 11))
+    expected_Hy[:-1, 1:, :-1] = np.tile(np.array([[[3]]], dtype=np.float32), (11, 11, 11))
+    expected_Hz[:-1, :-1, 1:] = np.tile(np.array([[[4]]], dtype=np.float32), (11, 11, 11))
+
+    # Why does fields_updates_normal use i+1, j+1 and k+1 everywhere?
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_magnetic()
+
+    assert np.equal(grid.Ex, expected_Ex).all()
+    assert np.equal(grid.Ey, expected_Ey).all()
+    assert np.equal(grid.Ez, expected_Ez).all()
+    assert np.equal(grid.Hx, expected_Hx).all()
+    assert np.equal(grid.Hy, expected_Hy).all()
+    assert np.equal(grid.Hz, expected_Hz).all()
+
+
+def test_update_electric_a_non_dispersive_zero_grid(config_mock):
+    grid = build_grid(100, 100, 100)
+
+    expected_value = np.zeros((101, 101, 101))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_electric_a()
+
+    assert np.equal(grid.Ex, expected_value).all()
+    assert np.equal(grid.Ey, expected_value).all()
+    assert np.equal(grid.Ez, expected_value).all()
+    assert np.equal(grid.Hx, expected_value).all()
+    assert np.equal(grid.Hy, expected_value).all()
+    assert np.equal(grid.Hz, expected_value).all()
+
+
+def test_update_electric_a_non_dispersive(config_mock):
+    grid = build_grid(11, 11, 11)
+
+    print(grid.updatecoeffsE)
+    print(grid.updatecoeffsE[1])
+    grid.updatecoeffsE[1] = 1
+    print(grid.updatecoeffsE[1])
+
+    grid.Ex = np.tile(np.array([[[3, 1], [1, 3]], [[1, 3], [3, 1]]], dtype=np.float32), (6, 6, 6))
+    grid.Ey = np.tile(np.array([[[1, 5], [5, 1]], [[5, 1], [1, 5]]], dtype=np.float32), (6, 6, 6))
+    grid.Ez = np.tile(np.array([[[5, 3], [3, 5]], [[3, 5], [5, 3]]], dtype=np.float32), (6, 6, 6))
+    grid.Hx = np.tile(np.array([[[1, 2], [2, 1]], [[2, 1], [1, 2]]], dtype=np.float32), (6, 6, 6))
+    grid.Hy = np.tile(np.array([[[1, 3], [3, 1]], [[3, 1], [1, 3]]], dtype=np.float32), (6, 6, 6))
+    grid.Hz = np.tile(np.array([[[1, 4], [4, 1]], [[4, 1], [1, 4]]], dtype=np.float32), (6, 6, 6))
+
+    expected_Ex = grid.Ex.copy()
+    expected_Ey = grid.Ey.copy()
+    expected_Ez = grid.Ez.copy()
+    expected_Hx = grid.Hx.copy()
+    expected_Hy = grid.Hy.copy()
+    expected_Hz = grid.Hz.copy()
+    # Why is there not a full (11x11x11) section of the frid being updated?
+    expected_Ex[:-1, 1:-1, 1:-1] = np.tile(np.array([[[2]]], dtype=np.float32), (11, 10, 10))
+    expected_Ey[1:-1, :-1, 1:-1] = np.tile(np.array([[[3]]], dtype=np.float32), (10, 11, 10))
+    expected_Ez[1:-1, 1:-1, :-1] = np.tile(np.array([[[4]]], dtype=np.float32), (10, 10, 11))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_electric_a()
+
+    assert np.equal(grid.Ex, expected_Ex).all()
+    assert np.equal(grid.Ey, expected_Ey).all()
+    assert np.equal(grid.Ez, expected_Ez).all()
+    assert np.equal(grid.Hx, expected_Hx).all()
+    assert np.equal(grid.Hy, expected_Hy).all()
+    assert np.equal(grid.Hz, expected_Hz).all()
+
+
+def test_update_electric_b_non_dispersive(config_mock):
+    grid = build_grid(100, 100, 100)
+
+    expected_value = np.zeros((101, 101, 101))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_electric_b()
+
+    assert np.equal(grid.Ex, expected_value).all()
+    assert np.equal(grid.Ey, expected_value).all()
+    assert np.equal(grid.Ez, expected_value).all()
+    assert np.equal(grid.Hx, expected_value).all()
+    assert np.equal(grid.Hy, expected_value).all()
+    assert np.equal(grid.Hz, expected_value).all()
+
+
+def test_update_electric_a_dispersive(config_mock):
+    assert False
+
+
+def test_update_electric_b_dispersive(config_mock):
+    assert False
+
+
+def test_update_magnetic_pml(config_mock):
+    grid = build_grid(100, 100, 100)
+
+    grid_expected_value = np.zeros((101, 101, 101))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_magnetic_pml()
+
+    assert np.equal(grid.Ex, grid_expected_value).all()
+    assert np.equal(grid.Ey, grid_expected_value).all()
+    assert np.equal(grid.Ez, grid_expected_value).all()
+    assert np.equal(grid.Hx, grid_expected_value).all()
+    assert np.equal(grid.Hy, grid_expected_value).all()
+    assert np.equal(grid.Hz, grid_expected_value).all()
+
+    for pml in grid.pmls["slabs"]:
+        assert np.equal(pml.HPhi1, 0).all()
+        assert np.equal(pml.HPhi2, 0).all()
+        assert np.equal(pml.EPhi1, 0).all()
+        assert np.equal(pml.EPhi2, 0).all()
+
+
+def test_update_electric_pml(config_mock):
+    grid = build_grid(100, 100, 100)
+
+    grid_expected_value = np.zeros((101, 101, 101))
+
+    cpu_updates = CPUUpdates(grid)
+    cpu_updates.update_electric_pml()
+
+    assert np.equal(grid.Ex, grid_expected_value).all()
+    assert np.equal(grid.Ey, grid_expected_value).all()
+    assert np.equal(grid.Ez, grid_expected_value).all()
+    assert np.equal(grid.Hx, grid_expected_value).all()
+    assert np.equal(grid.Hy, grid_expected_value).all()
+    assert np.equal(grid.Hz, grid_expected_value).all()
+
+    for pml in grid.pmls["slabs"]:
+        assert np.equal(pml.HPhi1, 0).all()
+        assert np.equal(pml.HPhi2, 0).all()
+        assert np.equal(pml.EPhi1, 0).all()
+        assert np.equal(pml.EPhi2, 0).all()
+
+
+def test_update_magnetic_sources(config_mock):
+    assert False
+
+
+def test_update_electric_sources(config_mock):
+    assert False
+
+
+def test_dispersive_update_a(config_mock):
+    assert False
+
+
+def test_dispersive_update_b(config_mock):
+    assert False

toolboxes/AntennaPatterns/initial_save.py

@@ -13,10 +13,15 @@ import h5py
 import matplotlib.pyplot as plt
 import numpy as np
 
-import gprMax.config as config
+from scipy.constants import c
+from scipy.constants import epsilon_0 as e0
+from scipy.constants import mu_0 as m0
+
 
 logger = logging.getLogger(__name__)
 
+# Impedance of free space (Ohms)
+z0 = np.sqrt(m0 / e0) 
 
 # Parse command line arguments
 parser = argparse.ArgumentParser(
@@ -59,9 +64,9 @@ traceno = np.s_[:]  # All traces
 # Critical angle and velocity
 if epsr:
     mr = 1
-    z1 = np.sqrt(mr / epsr) * config.sim_config.em_consts["z0"]
+    z1 = np.sqrt(mr / epsr) * z0
-    v1 = config.sim_config.em_consts["c"] / np.sqrt(epsr)
+    v1 = c / np.sqrt(epsr)
-    thetac = np.round(np.arcsin(v1 / config.sim_config.em_consts["c"]) * (180 / np.pi))
+    thetac = np.round(np.arcsin(v1 / c * (180 / np.pi))
     wavelength = v1 / f
 
 # Print some useful information
@@ -189,8 +194,8 @@ for radius in range(0, len(radii)):
             Ethetasum[index] = np.sum(Etheta[:, index] ** 2) / z1
             Hthetasum[index] = np.sum(Htheta[:, index] ** 2) / z1
         else:
-            Ethetasum[index] = np.sum(Etheta[:, index] ** 2) / config.sim_config.em_consts["z0"]
+            Ethetasum[index] = np.sum(Etheta[:, index] ** 2) / z0
-            Hthetasum[index] = np.sum(Htheta[:, index] ** 2) / config.sim_config.em_consts["z0"]
+            Hthetasum[index] = np.sum(Htheta[:, index] ** 2) / z0
 
         index += 1
 

toolboxes/AntennaPatterns/plot_fields.py

@@ -12,11 +12,17 @@ import os
 import matplotlib.pyplot as plt
 import numpy as np
 
-import gprMax.config as config
+from scipy.constants import c
+from scipy.constants import epsilon_0 as e0
+from scipy.constants import mu_0 as m0
+
 
 logger = logging.getLogger(__name__)
 
 
+# Impedance of free space (Ohms)
+z0 = np.sqrt(m0 / e0)
+
 # Parse command line arguments
 parser = argparse.ArgumentParser(
     description="Plot field patterns from a simulation with receivers positioned in circles around an antenna. This module should be used after the field pattern data has been processed and stored using the initial_save.py module.",
@@ -56,9 +62,9 @@ step = 12
 # Critical angle and velocity
 if epsr:
     mr = 1
-    z1 = np.sqrt(mr / epsr) * config.sim_config.em_consts["z0"]
+    z1 = np.sqrt(mr / epsr) * z0
-    v1 = config.sim_config.em_consts["c"] / np.sqrt(epsr)
+    v1 = c / np.sqrt(epsr)
-    thetac = np.round(np.rad2deg(np.arcsin(v1 / config.sim_config.em_consts["c"])))
+    thetac = np.round(np.rad2deg(np.arcsin(v1 / c)))
     wavelength = v1 / f
 
 # Print some useful information