Update examples level-set 2D and 3D.

Examples/levelSet2D.cl

@@ -5,18 +5,17 @@ __kernel void initScalar(__global float* scalar,
   uint gidX = get_global_id(0);
   uint gidY = get_global_id(1);
   uint i;
-  float r;
-  float2 pos, center=(float2)(0.5f, 0.75f);
+  float2 pos, center=(float2)(0.5, 0.75);
   for(i=gidX; i<WIDTH; i+=WGN)
     {
       pos = (float2)(i*size.x, gidY*size.y);
-      scalar[i+gidY*(WIDTH)] = ((distance(pos, center)<0.15f) ? 1.0f : 0.0f);
+      scalar[i+gidY*(WIDTH)] = ((distance(pos, center)<0.15) ? 1.0 : 0.0);
     }
 }
 
 
-__kernel void velocity(__global float* veloX,__global float* veloY,
-		       float t,
+__kernel void initVelocity(__global float* veloX,__global float* veloY,
+		       //float t,
 		       float4 minPos,
 		       float4 size)
 {
@@ -24,10 +23,10 @@ __kernel void velocity(__global float* veloX,__global float* veloY,
   uint gidY = get_global_id(1);
   uint i;
   float v;
-  float time_term = cos(t*M_PI_F/3.0f);
+  float time_term = 1.0;//cos(t*M_PI/3.0f);
   for(i=gidX; i<WIDTH; i+=WGN)
     {
-      v = veloX[i+gidY*(WIDTH)] = -sin(i*size.x * M_PI_F)*sin(i*size.x * M_PI_F)*sin(gidY*size.y * M_PI_F * 2)*time_term;
+      v = veloX[i+gidY*(WIDTH)] = -sin(i*size.x * M_PI)*sin(i*size.x * M_PI)*sin(gidY*size.y * M_PI * 2)*time_term;
       veloX[i+gidY*(WIDTH)] = v; // = Vx(x,y)
       veloY[i+gidY*(WIDTH)] = -v;// = Vy(y,x) = -Vx(x,y)
     }

Examples/levelSet2D.py

@@ -2,6 +2,7 @@
 import time
 import parmepy as pp
 import math
+import sys
 
 
 def vitesse(x, y):
@@ -18,18 +19,17 @@ def scalaire(x, y):
         return 0.
 
 
-def run():
+def run(nb=257):
     dim = 2
-    nb = 256
     boxLength = [1., 1.]
     boxMin = [0., 0.]
     nbElem = [nb, nb]
 
     timeStep = 0.05
     #period = 10.
-    finalTime = 3.
-    outputFilePrefix = './res/levelSet_2D_'
-    outputModulo = 1
+    finalTime = 40*timeStep
+    outputFilePrefix = './res_2D/levelSet_2D_'
+    outputModulo = 0
 
     t0 = time.time()
 
@@ -44,24 +44,27 @@ def run():
 
     ## Operators
     advec = pp.Transport(velo, scal)
-    velocity = pp.Velocity(velo)
+    #velocity = pp.Velocity(velo)
 
     ## Solver creation (discretisation of objects is done in solver initialisation)
     topo3D = pp.CartesianTopology(domain=box, resolution=nbElem, dim=dim)
 
     ##Problem
-    pb = pp.Problem(topo3D, [velocity, advec])
+    #pb = pp.Problem(topo3D, [velocity, advec])
+    pb = pp.Problem(topo3D, [advec])
 
     ## Setting solver to Problem
     pb.setSolver(finalTime, timeStep, 'gpu',
                  src=['./levelSet2D.cl'],
-                 io=pp.Printer(fields=[scal, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
+                 io=pp.Printer(fields=[scal, velo],
+                               frequency=outputModulo,
+                               outputPrefix=outputFilePrefix))
     pb.initSolver()
 
     t1 = time.time()
 
     ## Solve problem
-    pb.solve()
+    timings = pb.solve()
 
     tf = time.time()
 
@@ -69,6 +72,22 @@ def run():
     print "Total time : ", tf - t0, "sec (CPU)"
     print "Init time : ", t1 - t0, "sec (CPU)"
     print "Solving time : ", tf - t1, "sec (CPU)"
+    print ""
+    return pb.timings_info[0] + "\"Solving time\" \"Init time\" \"Total time\"", pb.timings_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0)
 
 if __name__ == "__main__":
-    run()
+    run(4097)
+    # timings = {}
+    # f = open('bench_levelSet2D.dat', 'w')
+    # header = "#size dim nPart "
+    # sizes = xrange(128, 4096+2048+1, 128)
+    # sizes = xrange(128, 512+1, 128)
+    # sizes = [1024,2048,4096,4096+2048]
+    # for s in sizes:
+    #     h, t = run(s)
+    #     timings[s] = t
+    # header += h
+    # f.write(header + "\n")
+    # for s in sizes:
+    #     f.write(str(s) + " " + str(2) + " " + str(s**2) + " " + str(timings[s]) + "\n")
+    # f.close()

Examples/levelSet3D.cl

@@ -7,37 +7,37 @@ __kernel void initScalar(__global float* scalar,
   uint gidZ = get_global_id(2);
   uint i;
   float r;
-  float4 pos, center=(float4)(0.35f, 0.35f, 0.35f, 0.0f);
+  float4 pos, center=(float4)(0.35, 0.35, 0.35, 0.0);
   for(i=gidX; i<WIDTH; i+=WGN)
     {
-      pos = (float4)(i*size.x, gidY*size.y, gidZ*size.z, 0.0f);
-      scalar[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = ((distance(pos, center)<0.15f) ? 1.0f : 0.0f);
+      pos = (float4)(i*size.x, gidY*size.y, gidZ*size.z, 0.0);
+      scalar[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = ((distance(pos, center)<0.15) ? 1.0 : 0.0);
     }
 }
 
-__kernel void velocity(__global float* veloX,
-		       __global float* veloY,
-		       __global float* veloZ,
-		       float t,
-		       float4 minPos,
-		       float4 size)
+__kernel void initVelocity(__global float* veloX,
+			   __global float* veloY,
+			   __global float* veloZ,
+			   //float t,
+			   float4 minPos,
+			   float4 size)
 {
   uint gidX = get_global_id(0);
   uint gidY = get_global_id(1);
   uint gidZ = get_global_id(2);
   uint i;
   float pix,piy,piz,v;
-  float time_term = cos(t*M_PI_F/3.0f);
+  float time_term = 1.0;//cos(t*M_PI/3.0f);
 
   for(i=gidX; i<WIDTH; i+=WGN)
     {
-      pix = i*size.x*M_PI_F;
-      piy = gidY*size.y*M_PI_F;
-      piz = gidZ*size.z*M_PI_F;
+      pix = i*size.x*M_PI;
+      piy = gidY*size.y*M_PI;
+      piz = gidZ*size.z*M_PI;
 
-      veloX[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = 2.0f * sin(pix)*sin(pix)*sin(2.0f*piy)*sin(2.0f*piz)*time_term; // Vx(x,y,z)
-      veloY[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0f*piy)*sin(pix)*sin(pix)*sin(2.0f*piz)*time_term;// Vy(y,x,z)
-      veloZ[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0f*piy)*sin(2.0f*piz)*sin(pix)*sin(pix)*time_term;//Vz(z,x,y)
+      veloX[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = 2.0 * sin(pix)*sin(pix)*sin(2.0*piy)*sin(2.0*piz)*time_term; // Vx(x,y,z)
+      veloY[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0*piy)*sin(pix)*sin(pix)*sin(2.0*piz)*time_term;// Vy(y,x,z)
+      veloZ[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0*piy)*sin(2.0*piz)*sin(pix)*sin(pix)*time_term;//Vz(z,x,y)
     }
 }
 

Examples/levelSet3D.py

@@ -2,20 +2,20 @@
 import time
 import parmepy as pp
 import math
+import sys
 
 
-def run():
+def run(nb=257):
     dim = 3
-    nb = 256
     boxLength = [1., 1., 1.]
     boxMin = [0., 0., 0.]
     nbElem = [nb, nb, nb]
 
     timeStep = 0.05
     #period = 10.
-    finalTime = 3.
-    outputFilePrefix = './res/levelSet_3D_'
-    outputModulo = 1
+    finalTime = 40*timeStep
+    outputFilePrefix = './res_3D/levelSet_3D_'
+    outputModulo = 0
 
     t0 = time.time()
 
@@ -30,13 +30,14 @@ def run():
 
     ## Operators
     advec = pp.Transport(velo, scal)
-    velocity = pp.Velocity(velo)
+    #velocity = pp.Velocity(velo)
 
     ## Solver creation (discretisation of objects is done in solver initialisation)
     topo3D = pp.CartesianTopology(domain=box, resolution=nbElem, dim=dim)
 
     ##Problem
-    pb = pp.Problem(topo3D, [velocity, advec])
+    #pb = pp.Problem(topo3D, [velocity, advec])
+    pb = pp.Problem(topo3D, [advec])
 
     ## Setting solver to Problem
     pb.setSolver(finalTime, timeStep, 'gpu',
@@ -47,7 +48,7 @@ def run():
     t1 = time.time()
 
     ## Solve problem
-    pb.solve()
+    timings = pb.solve()
 
     tf = time.time()
 
@@ -55,6 +56,21 @@ def run():
     print "Total time : ", tf - t0, "sec (CPU)"
     print "Init time : ", t1 - t0, "sec (CPU)"
     print "Solving time : ", tf - t1, "sec (CPU)"
+    print ""
+    return pb.timings_info[0] + "\"Solving time\" \"Init time\" \"Total time\"", pb.timings_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0)
 
 if __name__ == "__main__":
-    run()
+    run(257)
+    # timings = {}
+    # f = open('bench_levelSet3D.dat', 'w')
+    # header = "#size dim nPart "
+    # sizes = xrange(32, 256+32*3+1, 32)
+    # #sizes = [256+32+32+32]
+    # for s in sizes:
+    #     h, t = run(s)
+    #     timings[s] = t
+    # header += h
+    # f.write(header + "\n")
+    # for s in sizes:
+    #     f.write(str(s) + " " + str(3) + " " + str(s**3) + " " + str(timings[s]) + "\n")
+    # f.close()