start 1d transport test

test_dtw.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from __future__ import absolute_import, print_function
+import pyopencl as cl
+import numpy as np
+
+with open('kernels.cl', 'r') as f:
+    source = f.read()
+
+m = 2
+raf = 2
+deg = 2
+
+source = source.replace("THE_DEG", str(deg))
+source = source.replace("THE_RAF", str(raf))
+source = source.replace("THE_M", str(m))
+
+ctx = cl.create_some_context()
+
+prg = cl.Program(ctx, source).build()
+
+point_cpu = np.array([0.0, 0.0, 0.0,
+                  1.0, 0.0, 0.0,
+                  1.0, 1.0, 0.0,
+                  0.0, 1.0, 0.0,
+                  0.0, 0.0, 1.0,
+                  1.0, 0.0, 1.0,
+                  1.0, 1.0, 1.0,
+                  0.0, 1.0, 1.0], dtype = np.float32)
+
+
+
+element = [0,1,2,3,4,5,6,7]
+
+face2node = [[0,1,5,4],
+	      [1,2,6,5],
+	      [2,3,7,6],
+	      [0,4,7,3],
+	      [5,6,7,4],
+	      [0,3,2,1]]
+
+queue = cl.CommandQueue(ctx)
+mf = cl.mem_flags
+
+wsize = m * raf**3 * (deg+1)**3
+
+row_size = raf * (deg + 1)
+
+point_gpu = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf = point_cpu)
+wn_cpu = np.zeros(wsize, dtype = np.float32)
+wn_gpu = cl.Buffer(ctx, mf.WRITE_ONLY, size=(wsize * np.dtype('float32').itemsize))
+
+dtwn_cpu = np.zeros(wsize, dtype = np.float32)
+dtwn_gpu = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf = dtwn_cpu);
+
+ie = np.int32(0)
+
+prg.Init(queue, (row_size, row_size, row_size), None, ie, point_gpu, wn_gpu).wait()
+
+cl.enqueue_copy(queue, wn_cpu, wn_gpu)
+
+#print(wn_cpu)
+
+wn_out = wn_cpu.reshape((2*row_size**2,row_size))
+
+#print(wn_out)
+
+import matplotlib.pyplot as plt
+
+x_cpu = np.zeros(3 * row_size**3, dtype = np.float32)
+x_gpu = cl.Buffer(ctx, mf.WRITE_ONLY, x_cpu.nbytes)
+
+prg.get_nodes(queue, (row_size, row_size, row_size),
+              None, ie, point_gpu, x_gpu).wait()
+cl.enqueue_copy(queue, x_cpu, x_gpu).wait()
+
+plt.plot(x_cpu[0:row_size],wn_out[0])
+plt.plot(x_cpu[0:row_size],wn_out[1])
+
+#plt.show()
+#print(x_cpu[0:row_size])
+#print(x_cpu[0:row_size])
+
+tnow = np.float32(0)
+
+npts = deg+1
+
+prg.DGFlux(queue, (raf+1, row_size, row_size), (1, 1, 1), ie, point_gpu,tnow, wn_gpu,dtwn_gpu).wait()
+
+prg.DGVolume(queue, (row_size, row_size, row_size), (npts,npts,npts), ie, point_gpu, wn_gpu,dtwn_gpu).wait()
+
+cl.enqueue_copy(queue, dtwn_cpu, dtwn_gpu)
+
+print(dtwn_cpu)

testtranspose.py

@@ -3,6 +3,7 @@
 
 from __future__ import absolute_import, print_function
 import pyopencl as cl
+from pyopencl import array
 import numpy as np
 
 with open('kernels.cl', 'r') as f:
@@ -10,7 +11,7 @@ with open('kernels.cl', 'r') as f:
 
 m = 2
 raf = 2
-deg = 1
+deg = 2
 
 source = source.replace("THE_DEG", str(deg))
 source = source.replace("THE_RAF", str(raf))
@@ -48,23 +49,17 @@ wsize = m * raf**3 * (deg+1)**3
 row_size = raf * (deg + 1)
 
 point_gpu = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf = point_cpu)
-wn_cpu = np.zeros(wsize, dtype = np.float32)
-wn_gpu = cl.Buffer(ctx, mf.WRITE_ONLY, size=(wsize * np.dtype('float32').itemsize))
 
-dtwn_cpu = np.zeros(wsize, dtype = np.float32)
-dtwn_gpu = cl.Buffer(ctx, mf.READ_WRITE | mf.COPY_HOST_PTR, hostbuf = dtwn_cpu);
+wn = array.zeros(queue, shape=(wsize,), dtype=np.float32, order='C', allocator=None)
+wnp1 = array.zeros(queue, shape=(wsize,), dtype=np.float32, order='C', allocator=None)
 
-ie = np.int32(0)
-
-prg.Init(queue, (row_size, row_size, row_size), None, ie, point_gpu, wn_gpu).wait()
+dtwn = array.zeros(queue, shape=(wsize,), dtype=np.float32, order='C', allocator=None)
 
-cl.enqueue_copy(queue, wn_cpu, wn_gpu)
-
-#print(wn_cpu)
+ie = np.int32(0)
 
-wn_out = wn_cpu.reshape((2*row_size**2,row_size))
+prg.Init(queue, (row_size, row_size, row_size), None, ie, point_gpu, wn.base_data).wait()
 
-#print(wn_out)
+wn_out = wn.reshape((2*row_size**2,row_size))
 
 import matplotlib.pyplot as plt
 
@@ -73,10 +68,11 @@ x_gpu = cl.Buffer(ctx, mf.WRITE_ONLY, x_cpu.nbytes)
 
 prg.get_nodes(queue, (row_size, row_size, row_size),
               None, ie, point_gpu, x_gpu).wait()
+
 cl.enqueue_copy(queue, x_cpu, x_gpu).wait()
 
-plt.plot(x_cpu[0:row_size],wn_out[0])
-plt.plot(x_cpu[0:row_size],wn_out[1])
+#plt.plot(x_cpu[0:row_size],wn_out[0])
+#plt.plot(x_cpu[0:row_size],wn_out[1])
 
 #plt.show()
 #print(x_cpu[0:row_size])
@@ -86,10 +82,37 @@ tnow = np.float32(0)
 
 npts = deg+1
 
-prg.DGFlux(queue, (raf+1, row_size, row_size), (1, 1, 1), ie, point_gpu,tnow, wn_gpu,dtwn_gpu).wait()
+dx = 1./raf
+
+cfl = 0.1
+
+dt = dx * cfl
+
+t = 0.
+tmax = 0.05
+
+while t < tmax:
+
+    prg.DGFlux(queue, (raf+1, row_size, row_size), (1, 1, 1),
+               ie, point_gpu,tnow, wn.base_data,dtwn.base_data).wait()
+
+    prg.DGVolume(queue, (row_size, row_size, row_size), (npts,npts,npts),
+                 ie, point_gpu, wn.base_data,dtwn.base_data).wait()
+
+    wnp1 = wn + dt/2 * dtwn
+
+    dtwn.fill(0.)
+    
+    prg.DGFlux(queue, (raf+1, row_size, row_size), (1, 1, 1),
+               ie, point_gpu,tnow, wnp1.base_data,dtwn.base_data).wait()
+
+    prg.DGVolume(queue, (row_size, row_size, row_size), (npts,npts,npts),
+                 ie, point_gpu, wnp1.base_data,dtwn.base_data).wait()
+
+    wn = wn + dt * dtwn
 
-prg.DGVolume(queue, (row_size, row_size, row_size), (npts,npts,npts), ie, point_gpu, wn_gpu,dtwn_gpu).wait()
+    t = t + dt
+    print("t=",t)
 
-cl.enqueue_copy(queue, dtwn_cpu, dtwn_gpu)
 
-print(dtwn_cpu)
+print(dtwn)