Only initialise field/PML arrays on host when using CPU.

gprMax/grid.py

@@ -254,12 +254,12 @@ class FDTDGrid(Grid):
         import pycuda.gpuarray as gpuarray
 
         self.ID_gpu = gpuarray.to_gpu(self.ID)
-        self.Ex_gpu = gpuarray.to_gpu(self.Ex)
-        self.Ey_gpu = gpuarray.to_gpu(self.Ey)
-        self.Ez_gpu = gpuarray.to_gpu(self.Ez)
-        self.Hx_gpu = gpuarray.to_gpu(self.Hx)
-        self.Hy_gpu = gpuarray.to_gpu(self.Hy)
-        self.Hz_gpu = gpuarray.to_gpu(self.Hz)
+        self.Ex_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
+        self.Ey_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
+        self.Ez_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
+        self.Hx_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
+        self.Hy_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
+        self.Hz_gpu = gpuarray.to_gpu(np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=floattype))
 
     def gpu_initialise_dispersive_arrays(self):
         """Initialise dispersive material coefficient arrays on GPU."""

gprMax/model_build_run.py

@@ -180,7 +180,8 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
         G.initialise_geometry_arrays()
 
         # Initialise arrays for the field components
-        G.initialise_field_arrays()
+        if G.gpu is None:
+            G.initialise_field_arrays()
 
         # Process geometry commands in the order they were given
         process_geometrycmds(geometry, G)
@@ -296,12 +297,13 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
         if G.messages:
             print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
 
-        # Clear arrays for field components
-        G.initialise_field_arrays()
+        if G.gpu is None:
+            # Clear arrays for field components
+            G.initialise_field_arrays()
 
-        # Clear arrays for fields in PML
-        for pml in G.pmls:
-            pml.initialise_field_arrays()
+            # Clear arrays for fields in PML
+            for pml in G.pmls:
+                pml.initialise_field_arrays()
 
     # Adjust position of simple sources and receivers if required
     if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0:

gprMax/pml.py

@@ -196,7 +196,8 @@ class PML(object):
 
         self.CFS = G.cfs
 
-        self.initialise_field_arrays()
+        if G.gpu is None:
+            self.initialise_field_arrays()
 
     def initialise_field_arrays(self):
         """Initialise arrays to store fields in PML."""
@@ -308,19 +309,31 @@ class PML(object):
 
         import pycuda.gpuarray as gpuarray
 
-        self.EPhi1_gpu = gpuarray.to_gpu(self.EPhi1)
-        self.EPhi2_gpu = gpuarray.to_gpu(self.EPhi2)
         self.ERA_gpu = gpuarray.to_gpu(self.ERA)
         self.ERB_gpu = gpuarray.to_gpu(self.ERB)
         self.ERE_gpu = gpuarray.to_gpu(self.ERE)
         self.ERF_gpu = gpuarray.to_gpu(self.ERF)
-        self.HPhi1_gpu = gpuarray.to_gpu(self.HPhi1)
-        self.HPhi2_gpu = gpuarray.to_gpu(self.HPhi2)
         self.HRA_gpu = gpuarray.to_gpu(self.HRA)
         self.HRB_gpu = gpuarray.to_gpu(self.HRB)
         self.HRE_gpu = gpuarray.to_gpu(self.HRE)
         self.HRF_gpu = gpuarray.to_gpu(self.HRF)
 
+        if self.direction[0] == 'x':
+            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype))
+            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype))
+            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype))
+            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype))
+        elif self.direction[0] == 'y':
+            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype))
+            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=floattype))
+            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype))
+            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=floattype))
+        elif self.direction[0] == 'z':
+            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=floattype))
+            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=floattype))
+            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype))
+            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype))
+
     def gpu_get_update_funcs(self, kernelselectric, kernelsmagnetic):
         """Get update functions from PML kernels.
 
@@ -329,8 +342,6 @@ class PML(object):
             kernelsmagnetic: PyCuda SourceModule containing PML kernels for magnetic updates.
         """
 
-        from pycuda.compiler import SourceModule
-
         self.update_electric_gpu = kernelselectric.get_function('order' + str(len(self.CFS)) + '_' + self.direction)
         self.update_magnetic_gpu = kernelsmagnetic.get_function('order' + str(len(self.CFS)) + '_' + self.direction)