diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py
index b151cb9a..385f6f92 100644
--- a/gprMax/cmds_multiuse.py
+++ b/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)})"
diff --git a/gprMax/config.py b/gprMax/config.py
index 6d18bcad..d1d5385c 100644
--- a/gprMax/config.py
+++ b/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}"
- self.inputfilestr = Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
+ s = (
+ 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.scenes[model_num]
- except:
- return None
+ return sim_config.get_scene(self.model_num)
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.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)
+ self.general = {
+ "solver": "cpu",
+ "precision": "single",
+ "progressbars": args.log_level <= 20,
}
# 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 = args.scenes if args.scenes is not None else []
- except AttributeError:
- self.scenes = []
+ self.scenes: List[Optional[Scene]]
+ if hasattr(args, "scenes") and args.scenes is not None:
+ self.scenes = args.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):
- """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 get_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 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()
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index 7f94f79c..396dbad0 100644
--- a/gprMax/contexts.py
+++ b/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 .solvers import create_G, create_solver
+from .model import Model
+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
- can exploit parallelisation using either OpenMP (CPU), CUDA (GPU), or
- OpenCL (CPU/GPU).
+ """Standard context for building and running models.
+
+ 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.tsimstart = None
+ self.sim_start_time = 0
+ self.sim_end_time = 0
- def run(self):
- """Run the simulation in the correct context.
+ def _start_simulation(self) -> None:
+ """Run pre-simulation steps
- Returns:
- results: dict that can contain useful results/data from simulation.
+ Start simulation timer. Output copyright notice and host info.
"""
-
- self.tsimstart = timer()
+ self.sim_start_time = 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
- # called in a loop, and want to avoid this.
- config.model_configs = []
+ def _end_simulation(self) -> None:
+ """Run post-simulation steps
+
+ 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
- model_config = config.ModelConfig()
- config.model_configs.append(model_config)
+ self._run_model(i)
- # Always create a grid for the first model. The next model to run
- # 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()
+ self._end_simulation()
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):
- """Process for running a single model.
+ def _create_model_config(self, model_num: int) -> ModelConfig:
+ """Create model config and save to global config.
- Args:
- work: dict of any additional information that is passed to MPI
- workers. By default only model number (i) is used.
+ Set device in model config according to MPI rank.
"""
-
- # Create configuration for model
- config.model_num = work["i"]
- model_config = config.ModelConfig()
+ model_config = super()._create_model_config(model_num)
# 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()
- model = ModelBuildRun(G)
- model.build()
+ def _run_model(self, **work) -> None:
+ """Process for running a single model.
- if not config.sim_config.args.geometry_only:
- solver = create_solver(G)
- model.solve(solver)
- del solver, model
+ Args:
+ work: dict of any additional information that is passed to MPI
+ workers. By default only model number (i) is used.
+ """
+ return super()._run_model(work["i"])
- # Manual garbage collection required to stop memory leak on GPUs when
- # using pycuda
- del G
- gc.collect()
-
- def run(self):
+ def run(self) -> Optional[List[Optional[Dict]]]:
"""Specialise how the models are run.
Returns:
@@ -171,16 +238,12 @@ class TaskfarmContext(Context):
"""
if self.rank == 0:
- self.tsimstart = timer()
- 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"])
+ self._start_simulation()
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.print_sim_time_taken()
+ self._end_simulation()
return results
diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index 598e1dec..bc3fec28 100644
--- a/gprMax/gprMax.py
+++ b/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 .
-
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 "
- "can given in order to run 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.",
+ "scenes": (
+ "(list, req): Scenes to run the model. Multiple scene objects can given in order to run"
+ " 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.",
"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 "
- "typically be used to restart a series of models from a specific "
- "model number, with the n argument, e.g. to 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 "
- "models to be farmed out using a MPI task farm, e.g. to create a "
- "B-scan with 60 traces and use MPI to farm out each trace. For "
- "further details see the performance section of the User Guide.",
- "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).",
+ "i": (
+ "(int, opt): Model number to start/restart simulation from. It would typically be used to"
+ " restart a series of models from a specific model number, with the n argument, e.g. to"
+ " 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 models to be farmed out using a MPI task"
+ " farm, e.g. to create a B-scan with 60 traces and use MPI to farm out each trace. For"
+ " further details see the performance section of the User Guide."
+ ),
+ "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 "
- "with main grid coordinates to their equivalent local "
- "grid coordinate within the subgrid. If this option is "
- "off users must specify sub-grid object point within the "
- "global subgrid space.",
- "geometry_only": "(bool, opt): Build a model and produce any geometry " "views but do not run the simulation.",
- "geometry_fixed": "(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.",
+ "autotranslate": (
+ "(bool, opt): For sub-gridding - auto translate objects with main grid coordinates to their"
+ " equivalent local grid coordinate within the subgrid. If this option is off users must"
+ " specify sub-grid object point within the global subgrid space."
+ ),
+ "geometry_only": (
+ "(bool, opt): Build a model and produce any geometry views but do not run the simulation."
+ ),
+ "geometry_fixed": (
+ "(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
- simulation for the given list of scenes.
+ """Entry point for application programming interface (API).
+
+ Runs the simulation for the given list of scenes.
Args:
- scenes: list of the scenes to run the model. Multiple scene objects can
- be given in order to run multiple simulation runs. Each scene
- must contain the essential simulation objects.
- inputfile: optional string for input file path. Can also run simulation
- by providing an input file.
+ scenes: list of the scenes to run the model. Multiple scene
+ objects can be given in order to run multiple simulation
+ runs. Each scene must contain the essential simulation
+ objects.
+ 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.
- It would typically be used to restart a series of models from a
- specific model number, with the n argument, e.g. to restart from
- A-scan 45 when creating a B-scan with 60 traces.
- taskfarm: optional boolean flag to use Message Passing Interface (MPI) task
- farm. This option is most usefully combined with n to allow
- individual models to be farmed out using a MPI task farm,
- e.g. to create a B-scan with 60 traces and use MPI to farm out
- each trace. For further details see the performance section of
- the User Guide
- 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).
+ i: optional int for model number to start/restart simulation
+ from. It would typically be used to restart a series of
+ models from a specific model number, with the n argument,
+ e.g. to restart from A-scan 45 when creating a B-scan with
+ 60 traces.
+ taskfarm: optional boolean flag to use Message Passing Interface
+ (MPI) task farm. This option is most usefully combined with
+ n to allow individual models to be farmed out using a MPI
+ task farm, e.g. to create a B-scan with 60 traces and use
+ MPI to farm out each trace. For further details see the
+ performance section of the User Guide.
+ mpi: optional 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: 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
- objects with main grid coordinates to their equivalent
- local grid coordinate within the subgrid. If this option
- is off users must specify sub-grid object point within
- the global subgrid space.
+ autotranslate: optional boolean for sub-gridding to auto
+ translate objects with main grid coordinates to their
+ equivalent local grid coordinate within the subgrid. If this
+ option is off users must specify sub-grid object point
+ 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_fixed: optional boolean to run a series of models where the
- geometry does not change between models.
- write_processed: optional boolean to write another input file after any
- #python blocks (which are deprecated) in the
- original input file has been processed.
+ geometry views but do not run the simulation.
+ geometry_fixed: optional boolean to run a series of models where
+ the geometry does not change between models.
+ write_processed: optional boolean to write another input file
+ after any #python blocks (which are deprecated) in the
+ 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("--log-file", action="store_true", default=args_defaults["log_file"], help=help_msg["log_file"])
+ parser.add_argument(
+ "--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.
- Enables these to be propagated to calling script.
+ results: dict that can contain useful results/data from
+ 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()
diff --git a/gprMax/grid/cuda_grid.py b/gprMax/grid/cuda_grid.py
new file mode 100644
index 00000000..11adf31d
--- /dev/null
+++ b/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 .
+
+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)
diff --git a/gprMax/grid.py b/gprMax/grid/fdtd_grid.py
similarity index 52%
rename from gprMax/grid.py
rename to gprMax/grid/fdtd_grid.py
index 50693226..38b24640 100644
--- a/gprMax/grid.py
+++ b/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 .
-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 .utilities.utilities import fft_power, round_value
+# from gprMax.geometry_outputs import GeometryObjects, GeometryView
+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.hertziandipoles = []
- self.magneticdipoles = []
+ self.voltagesources: List[VoltageSource] = []
+ self.hertziandipoles: List[HertzianDipole] = []
+ self.magneticdipoles: List[MagneticDipole] = []
self.transmissionlines = []
- self.rxs = []
- self.srcsteps = [0, 0, 0]
- self.rxsteps = [0, 0, 0]
+ self.rxs: List[Rx] = []
+ self.srcsteps: List[float] = [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.Ey = np.zeros((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"])
+ self.Ex = np.zeros(
+ (self.nx + 1, self.ny + 1, self.nz + 1),
+ dtype=config.sim_config.dtypes["float_or_double"],
+ )
+ self.Ey = np.zeros(
+ (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.updatecoeffsH = np.zeros((len(self.materials), 5), dtype=config.sim_config.dtypes["float_or_double"])
+ self.updatecoeffsE = np.zeros(
+ (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"])))
diff --git a/gprMax/grid/mpi_grid.py b/gprMax/grid/mpi_grid.py
new file mode 100644
index 00000000..146aba93
--- /dev/null
+++ b/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 .
+
+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
diff --git a/gprMax/grid/opencl_grid.py b/gprMax/grid/opencl_grid.py
new file mode 100644
index 00000000..93a80e88
--- /dev/null
+++ b/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 .
+
+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)
diff --git a/gprMax/model.py b/gprMax/model.py
new file mode 100644
index 00000000..67adce00
--- /dev/null
+++ b/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 .
+
+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')}"
+ )
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
deleted file mode 100644
index b375eff7..00000000
--- a/gprMax/model_build_run.py
+++ /dev/null
@@ -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 .
-
-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)
diff --git a/gprMax/solvers.py b/gprMax/solvers.py
index 1e041468..91dab762 100644
--- a/gprMax/solvers.py
+++ b/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
-
-
-def create_G():
- """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
+from .updates.cpu_updates import CPUUpdates
+from .updates.cuda_updates import CUDAUpdates
+from .updates.opencl_updates import OpenCLUpdates
+from .updates.updates import Updates
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
diff --git a/gprMax/subgrids/grid.py b/gprMax/subgrids/grid.py
index 1711c2cb..1ac04fa3 100644
--- a/gprMax/subgrids/grid.py
+++ b/gprMax/subgrids/grid.py
@@ -17,8 +17,9 @@
# along with gprMax. If not, see .
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
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index 04805d93..a1cb4a67 100644
--- a/gprMax/subgrids/updates.py
+++ b/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
diff --git a/gprMax/taskfarm.py b/gprMax/taskfarm.py
index b9fde6b3..84f779d5 100644
--- a/gprMax/taskfarm.py
+++ b/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.")
diff --git a/gprMax/updates/cpu_updates.py b/gprMax/updates/cpu_updates.py
new file mode 100644
index 00000000..043e70a9
--- /dev/null
+++ b/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 .
+
+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
diff --git a/gprMax/updates/cuda_updates.py b/gprMax/updates/cuda_updates.py
new file mode 100644
index 00000000..c67b0d7a
--- /dev/null
+++ b/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 .
+
+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
diff --git a/gprMax/updates/opencl_updates.py b/gprMax/updates/opencl_updates.py
new file mode 100644
index 00000000..8438619b
--- /dev/null
+++ b/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 .
+
+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
diff --git a/gprMax/updates/updates.py b/gprMax/updates/updates.py
index 5da7e11a..182abb7f 100644
--- a/gprMax/updates/updates.py
+++ b/gprMax/updates/updates.py
@@ -1,1349 +1,96 @@
-# 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 .
-
-import logging
-from importlib import import_module
-
-import humanize
-import numpy as np
-from jinja2 import Environment, PackageLoader
-
-import gprMax.config as config
-
-from ..cuda_opencl import knl_fields_updates, knl_snapshots, knl_source_updates, knl_store_outputs
-from ..cython.fields_updates_normal import update_electric as update_electric_cpu
-from ..cython.fields_updates_normal import update_magnetic as update_magnetic_cpu
-from ..fields_outputs import store_outputs as store_outputs_cpu
-from ..receivers import dtoh_rx_array, htod_rx_arrays
-from ..snapshots import Snapshot, dtoh_snapshot_array, htod_snapshot_array
-from ..sources import htod_src_arrays
-from ..utilities.utilities import round32, timer
-
-logger = logging.getLogger(__name__)
-
-
-class CPUUpdates:
- """Defines update functions for CPU-based solver."""
-
- def __init__(self, G):
- """
- 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
-
-
-class CUDAUpdates:
- """Defines update functions for GPU-based (CUDA) solver."""
-
- def __init__(self, G):
- """
- 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
-
-
-class OpenCLUpdates:
- """Defines update functions for OpenCL-based solver."""
-
- def __init__(self, G):
- """
- 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
+# 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 .
+
+
+from abc import ABC, abstractmethod
+
+
+class Updates(ABC):
+ """Defines update functions for a solver."""
+
+ @abstractmethod
+ def store_outputs(self) -> None:
+ """Stores field component values for every receiver and transmission line."""
+ pass
+
+ @abstractmethod
+ def store_snapshots(self, iteration: int) -> None:
+ """Stores any snapshots.
+
+ Args:
+ iteration: int for iteration number.
+ """
+ pass
+
+ @abstractmethod
+ def update_magnetic(self) -> None:
+ """Updates magnetic field components."""
+ pass
+
+ @abstractmethod
+ def update_magnetic_pml(self) -> None:
+ """Updates magnetic field components with the PML correction."""
+ pass
+
+ @abstractmethod
+ def update_magnetic_sources(self) -> None:
+ """Updates magnetic field components from sources."""
+ pass
+
+ @abstractmethod
+ def update_electric_a(self) -> None:
+ """Updates electric field components."""
+ pass
+
+ @abstractmethod
+ def update_electric_pml(self) -> None:
+ """Updates electric field components with the PML correction."""
+ pass
+
+ @abstractmethod
+ def update_electric_sources(self) -> None:
+ """Updates electric field components from sources -
+ update any Hertzian dipole sources last.
+ """
+ pass
+
+ @abstractmethod
+ def update_electric_b(self) -> None:
+ """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.
+ """
+ pass
+
+ @abstractmethod
+ def time_start(self) -> None:
+ """Starts timer used to calculate solving time for model."""
+ pass
+
+ @abstractmethod
+ def calculate_solve_time(self) -> float:
+ """Calculates solving time for model."""
+ pass
+
+ def finalise(self) -> None:
+ pass
+
+ def cleanup(self) -> None:
+ pass
diff --git a/reframe_tests/job_scripts/archer2_tests.slurm b/reframe_tests/job_scripts/archer2_tests.slurm
index 63d36498..f3958fd2 100644
--- a/reframe_tests/job_scripts/archer2_tests.slurm
+++ b/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 "$@"
diff --git a/reframe_tests/regression_checks/BScanTest.h5 b/reframe_tests/regression_checks/BScanTest.h5
new file mode 100644
index 00000000..ff6c7602
Binary files /dev/null and b/reframe_tests/regression_checks/BScanTest.h5 differ
diff --git a/requirements.txt b/requirements.txt
index 646b4e9a..fcd838ba 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,7 +12,6 @@ matplotlib
numpy
numpy-stl
pandas
-pip
pre-commit
psutil
# pycuda
diff --git a/setup.py b/setup.py
index 743af3c6..409a97bd 100644
--- a/setup.py
+++ b/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(
diff --git a/tests/test_update_functions.py b/tests/test_update_functions.py
deleted file mode 100644
index 71541ab8..00000000
--- a/tests/test_update_functions.py
+++ /dev/null
@@ -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()
diff --git a/tests/updates/test_cpu_updates.py b/tests/updates/test_cpu_updates.py
new file mode 100644
index 00000000..73fb93fb
--- /dev/null
+++ b/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
diff --git a/toolboxes/AntennaPatterns/initial_save.py b/toolboxes/AntennaPatterns/initial_save.py
index 1c849ab4..7a7fe4cb 100644
--- a/toolboxes/AntennaPatterns/initial_save.py
+++ b/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"]
- v1 = config.sim_config.em_consts["c"] / np.sqrt(epsr)
- thetac = np.round(np.arcsin(v1 / config.sim_config.em_consts["c"]) * (180 / np.pi))
+ z1 = np.sqrt(mr / epsr) * z0
+ v1 = c / np.sqrt(epsr)
+ 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"]
- Hthetasum[index] = np.sum(Htheta[:, index] ** 2) / config.sim_config.em_consts["z0"]
+ Ethetasum[index] = np.sum(Etheta[:, index] ** 2) / z0
+ Hthetasum[index] = np.sum(Htheta[:, index] ** 2) / z0
index += 1
diff --git a/toolboxes/AntennaPatterns/plot_fields.py b/toolboxes/AntennaPatterns/plot_fields.py
index d6f85e4b..b43160ed 100644
--- a/toolboxes/AntennaPatterns/plot_fields.py
+++ b/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"]
- v1 = config.sim_config.em_consts["c"] / np.sqrt(epsr)
- thetac = np.round(np.rad2deg(np.arcsin(v1 / config.sim_config.em_consts["c"])))
+ z1 = np.sqrt(mr / epsr) * z0
+ v1 = c / np.sqrt(epsr)
+ thetac = np.round(np.rad2deg(np.arcsin(v1 / c)))
wavelength = v1 / f
# Print some useful information