Fix gpu solver for data transfers

Examples/mainJM.py

@@ -33,6 +33,8 @@ def run():
     ## Fields
     scal = pp.ContinuousField(domain=box, name='Scalar')
     velo = pp.ContinuousField(domain=box, name='Velocity', vector=True)
+    #scal = pp.AnalyticalField(domain=box, name='Scalar')
+    #velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
 
     ## Operators
     advec = pp.Transport(velo, scal)
@@ -44,10 +46,10 @@ def run():
     pb = pp.Problem(topo3D, [advec])
 
     ## Setting solver to Problem
-    pb.setSolver(finalTime, timeStep, 'gpu', src='./mainJM_kernels.cl')
-
-    # ## Setting io to Problem
-    pb.setIO(pp.Printer(fields=[scal, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
+    #pb.setSolver(finalTime, timeStep, 'gpu',
+    pb.setSolver(finalTime, timeStep, 'gpu', src='./mainJM_kernels.cl',
+                 io=pp.Printer(fields=[scal, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
+    pb.initSolver()
 
     t1 = time.time()
 

HySoP/hysop/fields/discrete.py

@@ -71,7 +71,7 @@ class ScalarField(object):
                 self.data = v_formula(self.topology.mesh.coords[0])
             self.contains_data = True
         else:
-            print "",
+            print "No formula",
 
     def tofile(self, filename):
         evtk.imageToVTK(filename, pointData={self.name: self.data})

HySoP/hysop/operator/advection_d.py

@@ -81,6 +81,9 @@ class Advection_P(DiscreteOperator):
                                         coord_min,
                                         domain_length,
                                         mesh_size)
+            # Output data of a OpenCL kernels are on device, host data are outdated
+            for df in self.output:
+                df.contains_data = False
             self.numMethod.queue.finish()
             c_time = (evt_copy.profile.end - evt_copy.profile.start) * 1e-9
             c_time += (evt.profile.end - evt.profile.start) * 1e-9

HySoP/hysop/operator/remeshing.py

@@ -81,6 +81,9 @@ class Remeshing(DiscreteOperator):
                                         PARMES_INTEGER_GPU(nb_ligne),
                                         coord_min,
                                         mesh_size)
+            # Output data of a OpenCL kernels are on device, host data are outdated
+            for df in self.output:
+                df.contains_data = False
             self.numMethod.queue.finish()
             c_time = (evt.profile.end - evt.profile.start) * 1e-9
         else:

HySoP/hysop/particular_solvers/basic.py

@@ -39,6 +39,7 @@ class ParticleSolver(Solver):
         self.RemeshingMethod = RemeshingMethod
         ## Splitting Method
         self.splittingMethod = splittingMethod
+        self.isInitialized = False
         if timer is None:
             self.problem.setTimer(Timer(t_end, dt))
         else:
@@ -90,6 +91,8 @@ class ParticleSolver(Solver):
                     self.problem.operators[index] = Splitting([op, self.remesh])
                 else:
                     self.problem.operators[index] = Splitting([op])
+        self.problem.io.init()
+        self.isInitialized = True
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/particular_solvers/gpu.py

@@ -13,6 +13,7 @@ from ..tools.timer import Timer
 from ..tools.printer import Printer
 import pyopencl as cl
 from ..tools.gpu_data_transfer import hostToDevice
+from ..tools.gpu_data_transfer import deviceToHost
 
 
 class GPUParticleSolver(ParticleSolver):
@@ -34,7 +35,13 @@ class GPUParticleSolver(ParticleSolver):
         Constructor.
         Create a solver description.
         """
-        ParticleSolver.__init__(self, problem, t_end, dt)
+        ParticleSolver.__init__(self, problem, t_end, dt,
+                                ODESolver=ODESolver,
+                                InterpolationMethod=InterpolationMethod,
+                                RemeshingMethod=RemeshingMethod,
+                                splittingMethod=splittingMethod,
+                                timer=timer,
+                                io=io)
         self.user_gpu_src = src
         print "=== Discover OpenCL environment ==="
         #Get platform.
@@ -49,7 +56,14 @@ class GPUParticleSolver(ParticleSolver):
         print "  - Name       :", self.platform.name
         print "  - Version    :", self.platform.version
         #Get device.
-        self.device = self.platform.get_devices(eval("cl.device_type." + device_type.upper()))[device_id]
+        try:
+            self.device = self.platform.get_devices(eval("cl.device_type." + device_type.upper()))[device_id]
+        except cl.RuntimeError as e:
+            print "RuntimeError:", e
+            self.device = cl.create_some_context().devices[0]
+        except AttributeError as e:
+            print "AttributeError:", e
+            self.device = cl.create_some_context().devices[0]
         print "  Device"
         print "  - Name                :", self.device.name
         print "  - Type                :", cl.device_type.to_string(self.device.type)
@@ -144,14 +158,14 @@ class GPUParticleSolver(ParticleSolver):
                             evt = initKernel(self.queue,
                                              global_wg,
                                              local_wg,
-                                             df.gpu_data[0], df.gpu_data[1], df.gpu_data[2],
+                                             df.gpu_data[XDIR], df.gpu_data[YDIR], df.gpu_data[ZDIR],
                                              coord_min,
                                              mesh_size)
                         else:
                             evt = initKernel(self.queue,
                                              global_wg,
                                              local_wg,
-                                             df.gpu_data[0], df.gpu_data[1],
+                                             df.gpu_data[XDIR], df.gpu_data[YDIR],
                                              coord_min,
                                              mesh_size)
                     else:
@@ -179,27 +193,11 @@ class GPUParticleSolver(ParticleSolver):
         ## Setting advection and remeshing kernels:
         self.advection.setMethod(KernelLauncher(self.prg.advection, self.queue))
         self.remesh.setMethod(KernelLauncher(self.prg.remeshing, self.queue))
-
-        # total_data, total_time = 0, 0.
-        # for f in self.problem.variables:
-        #     initKernel = None
-        #     for k in self.prg.all_kernels():
-        #         if k.get_info(cl.kernel_info.FUNCTION_NAME) == ('init'+f.name):
-        #             initKernel = k
-        #     if initKernel is not None:
-        #         print "Initialisation kernel found", initKernel
-        #     else:
-        #         for df in f.discreteField:
-        #             if df.contains_data:
-        #                 print f.name, ":",
-        #                 temp_data, temp_time = hostToDevice(self.queue, df)
-        #         total_data += temp_data
-        #         total_time += temp_time
-        # if total_data > 0:
-        #     print "Total : ", total_data, "Bytes transfered at {0:.3f} GBytes/sec".format((total_data * 1e-9) / total_time)
-        # self.queue.finish()
+        self.problem.io.set_get_data_method(self.get_data_from_device)
         print "===\n"
-        ##\todo brancher les objets opencl sur les opérateurs et les fields
+
+    def get_data_from_device(self, field):
+        deviceToHost(self.queue, field)
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/problem/problem.py

@@ -43,10 +43,16 @@ class Problem():
             self.solver = GPUParticleSolver(self, t_end, dt, **kwargs)
         else:
             raise ValueError("Unknown solver type : " + str(solver_type))
-        self.solver.initialize()
+
+    def initSolver(self):
+        if self.solver is not None:
+            self.solver.initialize()
 
     def solve(self):
         print "\n\n Start solving ..."
+        if not self.solver.isInitialized:
+            self.initSolver()
+        self.io.step()
         while not self.timer.end:
             print "==== Iteration : {0:3d}   t={1:6.2f} ====".format(self.timer.ite + 1, self.timer.t + self.timer.dt)
             for op in self.operators:

HySoP/hysop/tools/gpu_data_transfer.py

@@ -5,6 +5,7 @@ GPU Memory transfering methods
 """
 from ..constants import *
 import pyopencl as cl
+import numpy as np
 
 
 def hostToDevice(queue, discreteField):
@@ -34,24 +35,35 @@ def hostToDevice(queue, discreteField):
 
 def deviceToHost(queue, discreteField):
     print "device->host :",
-    data, time = 0, 0.
-    # if discreteField.vector:
-    #     evt = cl.enqueue_copy(queue, discreteField.gpu_data[XDIR],
-    #                           discreteField.data[XDIR].ravel())
-    #     data += discreteField.data[XDIR].nbytes
-    #     time += (evt.profile.end - evt.profile.start) * 1e-9
-    #     evt = cl.enqueue_copy(queue, discreteField.gpu_data[YDIR],
-    #                           discreteField.data[YDIR].swapaxes(1, 0).ravel())
-    #     data += discreteField.data[YDIR].nbytes
-    #     time += (evt.profile.end - evt.profile.start) * 1e-9
-    #     if len(discreteField.data) == 3:
-    #         evt = cl.enqueue_copy(queue, discreteField.gpu_data[ZDIR],
-    #                               discreteField.data[ZDIR].swapaxes(1, 0).swapaxes(2, 0).ravel())
-    #         data += discreteField.data[ZDIR].nbytes
-    #         time += (evt.profile.end - evt.profile.start) * 1e-9
-    # else:
-    #     evt = cl.enqueue_copy(queue, discreteField.gpu_data, discreteField.data)
-    #     data += discreteField.data.nbytes
-    #     time += (evt.profile.end - evt.profile.start) * 1e-9
-    print data, "Bytes transfered at {0:.3f} GBytes/sec".format(0)  # (data * 1e-9) / time)
+    data0, data1, data2 = 0, 0, 0
+    time0, time1, time2 = 0., 0., 0.
+    if discreteField.vector:
+        # Use of discreteField.data[XDIR] as temporary array \todo Need cleaner solution
+        # Get ZDIR component
+        if len(discreteField.data) == 3:
+            evt = cl.enqueue_copy(queue, discreteField.data[XDIR],
+                                  discreteField.gpu_data[ZDIR])
+            discreteField.data[ZDIR] = discreteField.data[XDIR].swapaxes(2, 0).swapaxes(1, 0)
+            data0 = discreteField.data[ZDIR].nbytes
+            time0 = (evt.profile.end - evt.profile.start) * 1e-9
+        # Get YDIR component
+        evt = cl.enqueue_copy(queue, discreteField.data[XDIR],
+                              discreteField.gpu_data[YDIR])
+        discreteField.data[YDIR] = discreteField.data[XDIR].swapaxes(1, 0)
+        data1 = discreteField.data[YDIR].nbytes
+        time1 = (evt.profile.end - evt.profile.start) * 1e-9
+        # Get XDIR component
+        evt = cl.enqueue_copy(queue, discreteField.data[XDIR],
+                              discreteField.gpu_data[XDIR])
+        data2 = discreteField.data[XDIR].nbytes
+        time2 = (evt.profile.end - evt.profile.start) * 1e-9
+    else:
+        evt = cl.enqueue_copy(queue, discreteField.data,
+                              discreteField.gpu_data)
+        data0 = discreteField.data.nbytes
+        time0 = (evt.profile.end - evt.profile.start) * 1e-9
+    queue.finish()
+    data = data0 + data1 + data2
+    time = time0 + time1 + time2
+    print data, "Bytes transfered at {0:.3f} GBytes/sec".format(data / (time * 1024 ** 3))
     return data, time

HySoP/hysop/tools/printer.py

@@ -25,10 +25,13 @@ class Printer():
         self.outputPrefix = outputPrefix
         self.fields = fields
         self.ite = 0
+        self.get_data_method = None
         if self.freq != 0:
             self.step = self._printStep
         else:
             self.step = self._passStep
+
+    def init(self):
         self.__dict_data = {}
         for f in self.fields:
             for df in f.discreteField:
@@ -37,17 +40,30 @@ class Printer():
                         self.__dict_data[df.name + S_DIR[d]] = df.data[d]
                 else:
                     self.__dict_data[df.name] = df.data
-        self.step()
+
+    def set_get_data_method(self, method):
+        if callable(method):
+            self.get_data_method = method
+        else:
+            raise ValueError("Cannot set non callable method to get data to print. Given method : " + str(method))
 
     def _passStep(self):
         pass
 
     def _printStep(self):
+        print "== IO"
+        if self.get_data_method is not None:
+            for f in self.fields:
+                for df in f.discreteField:
+                    if not df.contains_data:
+                        print "Data transfer : ", f.name, ":",
+                        self.get_data_method(df)
         if (self.ite % self.freq) == 0:
             filename = self.outputPrefix + "results_{0:05d}.dat".format(self.ite)
             print "Print file : " + filename
             evtk.imageToVTK(filename, pointData=self.__dict_data)
         self.ite += 1
+        print "==\n"
     #     self.skip -= 1
     #     if self.skip == 0:
     #         if isinstance(self.problem.discreteProblem.solver, GPUParticularSolver) and self.ite > 0: