diff --git a/HySoP/CMake/ParmesTests.cmake b/HySoP/CMake/ParmesTests.cmake
index 5f1bd553902eea89fa03a458db3bcae540c941ff..5888dad433261f14c4f2f71db1576bef6c5efeac 100644
--- a/HySoP/CMake/ParmesTests.cmake
+++ b/HySoP/CMake/ParmesTests.cmake
@@ -8,7 +8,7 @@ find_python_module(pytest REQUIRED)
set(py_src_dirs
fields
domain
- operator
+ #operator
problem
tools
numerics
@@ -19,6 +19,11 @@ if(USE_MPI)
${py_src_dirs} mpi)
endif()
+if(WITH_GPU)
+ set(py_src_dirs
+ ${py_src_dirs} gpu)
+endif()
+
### JM : Bub, GLOB_RECURSE utilse les chemins relatifs avec ../
# foreach(testdir ${py_src_dirs})
# file(GLOB_RECURSE py_test_files "${test_dir}" RELATIVE_PATH parmepy/*/test_*.py)
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
index 75e0b3b263842f5c97a151594c203dd4a118a6d0..9fc1a1d16f72b88e2e86b0613fb75015149c453f 100644
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -46,7 +46,7 @@ if __DEBUG__:
print ' (Inherits from : ',
print [c.__name__ for c in args[0].__mro__[1:]], ')'
else:
- print 'Call : ', f.__name__, 'of ', f.__module__
+ print 'Call : ', f.__name__, 'in ', f.__code__.co_filename, f.__code__.co_firstlineno
## Calling f
r = f(*args, **kw)
if f.__name__ is '__new__':
diff --git a/HySoP/hysop/domain/domain.py b/HySoP/hysop/domain/domain.py
index 7cfd811082e1fd82b1773fed3e1449204d6bebff..c1e85e5735a64b48568c41b09b1fd3d23035243b 100644
--- a/HySoP/hysop/domain/domain.py
+++ b/HySoP/hysop/domain/domain.py
@@ -23,11 +23,19 @@ class Domain(object):
"""
## Domain dimension.
self.dimension = dimension
- ## Domain topologies
+ ## Domain topologies. List of uniques topologies.
self.topologies = []
@debug
def addTopology(self, topo):
+ """Topology management.
+ If an equivalent topology is found in topologies, then return this
+ topology from the list and delete the other.
+ Otherwise, add the topology in the list.
+
+ @param topo : Topology to add
+ @return topology and its index in the topology list
+ """
try:
idTopo = self.topologies.index(topo)
del topo
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
index c47729ff3f26680c511f710f3f2a50c2e7297f24..acb8fd772b5618e365672f7c0a9460ebb39e2183 100644
--- a/HySoP/hysop/fields/analytical.py
+++ b/HySoP/hysop/fields/analytical.py
@@ -14,7 +14,7 @@ class AnalyticalField(Field):
"""
@debug
- def __init__(self, domain, formula, name="", vector=False):
+ def __init__(self, domain, formula, name="", isVector=False):
"""
Create an AnalyticalField from a formula.
@@ -23,9 +23,10 @@ class AnalyticalField(Field):
@param name : name of the variable (used for vtk output).
@param vector : is variable is a vector.
- @note formula is used in ScalarField or VectorField as vectorized function by numpy.
+ @note formula is used in ScalarField or VectorField as vectorized
+ function by numpy.
"""
- Field.__init__(self, domain, name, vector)
+ Field.__init__(self, domain, name, isVector)
## Analytic formula.
self.formula = formula
@@ -41,7 +42,8 @@ class AnalyticalField(Field):
@debug
def initialize(self):
if self._fieldId == -1:
- raise ValueError("Cannot initialise analytical field non discretized.")
+ raise ValueError(
+ "Cannot initialise analytical field non discretized.")
print self.name,
for dF in self.discreteField:
dF.initialize(self.formula)
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 61d83ca9b504fe90294b5d28bf27772eed617e78..6d1f8612e079620d8267dd5834772b3167b07de1 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -1,7 +1,7 @@
"""
@package parmepy.fields.continuous
-Continuous variable description
+Continuous variable description.
"""
from parmepy import __VERBOSE__
from parmepy.constants import debug
@@ -22,15 +22,13 @@ class Field(object):
Create a Field.
@param domain : variable application domain.
- @param name : name of the variable (used for vtk output).
+ @param name : name of the variable.
@param vector : is variable is a vector.
"""
## Application domain of the variable.
self.domain = domain
## Field dimension.
self.dimension = self.domain.dimension
- ## Discretisation of this field
- self.discreteField = None
## Name of this field
if(name is not None):
self.name = name
@@ -49,9 +47,11 @@ class Field(object):
def discretize(self, topo):
"""
Discretization of the field on a topology.
+ In the case of this field already discretized whith the same topology,
+ a discrete object is not created.
+
@param topo : topology specification for discretization
@return discrete field and its index.
-
"""
# Find if topology already use to discretize self
try:
@@ -74,7 +74,6 @@ class Field(object):
name=nameD,
idFromParent=fieldId))
self.topologies.insert(fieldId, topo)
-
return self.discreteField[fieldId], fieldId
@debug
diff --git a/HySoP/hysop/fields/gpu_discrete.py b/HySoP/hysop/fields/gpu_discrete.py
index c7fb28ca2aca7994b27d3d288a1e631a56aeeb90..0e2a7b023fa9d17e197a4642adfd2220b2b1c134 100644
--- a/HySoP/hysop/fields/gpu_discrete.py
+++ b/HySoP/hysop/fields/gpu_discrete.py
@@ -20,6 +20,7 @@ class GPUDiscreteField(object):
self.initialization_kernel = None
self.precision = precision
self.mem_size = 0
+ self.data_on_device = False
@abstractmethod
def toHost(self):
@@ -74,47 +75,48 @@ class GPUVectorField(VectorField, GPUDiscreteField):
@debug
def initialize(self, formula=None):
- if self.initialization_kernel is None:
- VectorField.initialize(self, formula)
- if self.dimension == 2:
- self.data[0] = np.asarray(self.data[0],
- dtype=self.precision, order=ORDER)
- self.data[1] = np.asarray(self.data[1],
- dtype=self.precision, order=ORDER)
- else:
- self.data[0] = np.asarray(self.data[0],
- dtype=self.precision, order=ORDER)
- self.data[1] = np.asarray(self.data[1],
- dtype=self.precision, order=ORDER)
- self.data[2] = np.asarray(self.data[2],
- dtype=self.precision, order=ORDER)
- data, time = self.toDevice()
- print "Transfers : ", self.mem_size, "Bytes transfered at ",
- print "{0:.3f} GBytes/sec".format(
- (self.mem_size * 1e-9) / time)
- else:
- coord_min = np.ones(4, dtype=self.precision)
- mesh_size = np.ones(4, dtype=self.precision)
- coord_min[0:self.dimension] = np.asarray(self.topology.mesh.origin,
- dtype=self.precision)
- mesh_size[0:self.dimension] = np.asarray(self.topology.mesh.space_step,
- dtype=self.precision)
- if self.dimension == 2:
- evt = self.initialization_kernel(self.gpu_data[0],
- self.gpu_data[1],
- coord_min,
- mesh_size)
+ if not self.data_on_device:
+ if self.initialization_kernel is None:
+ VectorField.initialize(self, formula)
+ if self.dimension == 2:
+ self.data[0] = np.asarray(self.data[0],
+ dtype=self.precision, order=ORDER)
+ self.data[1] = np.asarray(self.data[1],
+ dtype=self.precision, order=ORDER)
+ else:
+ self.data[0] = np.asarray(self.data[0],
+ dtype=self.precision, order=ORDER)
+ self.data[1] = np.asarray(self.data[1],
+ dtype=self.precision, order=ORDER)
+ self.data[2] = np.asarray(self.data[2],
+ dtype=self.precision, order=ORDER)
+ data, time = self.toDevice()
+ print "Transfers : ", self.mem_size, "Bytes transfered at ",
+ print "{0:.3f} GBytes/sec".format(
+ (self.mem_size * 1e-9) / time)
else:
- evt = self.initialization_kernel(self.gpu_data[0],
- self.gpu_data[1],
- self.gpu_data[2],
- coord_min,
- mesh_size)
- self.queue.finish()
- cTime = (evt.profile.end - evt.profile.start) * 1e-9
- print "Computing : ", cTime, "sec"
- self.contains_data = False
- self.data_on_device = True
+ coord_min = np.ones(4, dtype=self.precision)
+ mesh_size = np.ones(4, dtype=self.precision)
+ coord_min[0:self.dimension] = np.asarray(self.topology.mesh.origin,
+ dtype=self.precision)
+ mesh_size[0:self.dimension] = np.asarray(self.topology.mesh.space_step,
+ dtype=self.precision)
+ if self.dimension == 2:
+ evt = self.initialization_kernel(self.gpu_data[0],
+ self.gpu_data[1],
+ coord_min,
+ mesh_size)
+ else:
+ evt = self.initialization_kernel(self.gpu_data[0],
+ self.gpu_data[1],
+ self.gpu_data[2],
+ coord_min,
+ mesh_size)
+ self.queue.finish()
+ cTime = (evt.profile.end - evt.profile.start) * 1e-9
+ print "Computing : ", cTime, "sec"
+ self.contains_data = False
+ self.data_on_device = True
def toDevice(self):
"""
@@ -245,29 +247,30 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
@debug
def initialize(self, formula=None):
- if self.initialization_kernel is None:
- ScalarField.initialize(self, formula)
- self.data = np.asarray(self.data,
- dtype=self.precision, order=ORDER)
- data, time = self.toDevice()
- print "Transfers : ", self.mem_size, "Bytes transfered at ",
- print "{0:.3f} GBytes/sec".format(
- (self.mem_size * 1e-9) / time)
- else:
- coord_min = np.ones(4, dtype=self.precision)
- mesh_size = np.ones(4, dtype=self.precision)
- coord_min[0:self.dimension] = np.asarray(self.topology.mesh.origin,
- dtype=self.precision)
- mesh_size[0:self.dimension] = np.asarray(self.topology.mesh.space_step,
- dtype=self.precision)
- evt = self.initialization_kernel(self.gpu_data,
- coord_min,
- mesh_size)
- self.queue.finish()
- cTime = (evt.profile.end - evt.profile.start) * 1e-9
- print "Computing : ", cTime, "sec"
- self.contains_data = False
- self.data_on_device = True
+ if not self.data_on_device:
+ if self.initialization_kernel is None:
+ ScalarField.initialize(self, formula)
+ self.data = np.asarray(self.data,
+ dtype=self.precision, order=ORDER)
+ data, time = self.toDevice()
+ print "Transfers : ", self.mem_size, "Bytes transfered at ",
+ print "{0:.3f} GBytes/sec".format(
+ (self.mem_size * 1e-9) / time)
+ else:
+ coord_min = np.ones(4, dtype=self.precision)
+ mesh_size = np.ones(4, dtype=self.precision)
+ coord_min[0:self.dimension] = np.asarray(self.topology.mesh.origin,
+ dtype=self.precision)
+ mesh_size[0:self.dimension] = np.asarray(self.topology.mesh.space_step,
+ dtype=self.precision)
+ evt = self.initialization_kernel(self.gpu_data,
+ coord_min,
+ mesh_size)
+ self.queue.finish()
+ cTime = (evt.profile.end - evt.profile.start) * 1e-9
+ print "Computing : ", cTime, "sec"
+ self.contains_data = False
+ self.data_on_device = True
def toDevice(self):
"""
diff --git a/HySoP/hysop/fields/scalar.py b/HySoP/hysop/fields/scalar.py
index 2887f6f83883ce97e4dc22a8acd4568906cc8bab..9c811fe5b3b959320e804c6e35b1678649798255 100644
--- a/HySoP/hysop/fields/scalar.py
+++ b/HySoP/hysop/fields/scalar.py
@@ -57,7 +57,7 @@ class ScalarField(DiscreteField):
@li 2D : float formula(float x, float y)
where x,y,z are point coordinates.
"""
- if formula is not None:
+ if formula is not None and not self.contains_data:
if __VERBOSE__:
print "...",
v_formula = np.vectorize(formula)
diff --git a/HySoP/hysop/fields/vector.py b/HySoP/hysop/fields/vector.py
index 6451f34e85c89c54b6aa539f489068e2374d582c..313b3dbc07b39d375574575662d0eb8540602dd4 100644
--- a/HySoP/hysop/fields/vector.py
+++ b/HySoP/hysop/fields/vector.py
@@ -87,7 +87,7 @@ class VectorField(DiscreteField):
@li 2D : float, float formula(float x, float y)
where x,y,z are point coordinates.
"""
- if formula is not None:
+ if formula is not None and not self.contains_data:
if __VERBOSE__:
print "...",
v_formula = np.vectorize(formula)
diff --git a/HySoP/hysop/gpu/cl_src/weights_m4prime.cl b/HySoP/hysop/gpu/cl_src/weights_m4prime.cl
index 866122bb746dfb17273090124f8cf3f13feb45a1..0a3337eef58bde74723f6613b18aed61d7dd0f80 100644
--- a/HySoP/hysop/gpu/cl_src/weights_m4prime.cl
+++ b/HySoP/hysop/gpu/cl_src/weights_m4prime.cl
@@ -7,10 +7,11 @@
* @return weights
*/
inline float alpha(float y, float s){
- return ((y * (y * (-y + 2.0) - 1.0)) / 2.0) * s;}
+ return ((y * (y * (-y + 2.0) - 1.0)) / 2.0) * s;} //6Flops
inline float beta(float y, float s){
- return ((y * y * (3.0 * y - 5.0) + 2.0) / 2.0) * s;}
+ return ((y * y * (3.0 * y - 5.0) + 2.0) / 2.0) * s;} //7flops
inline float gamma(float y, float s){
- return ((y * (y * (-3.0 * y + 4.0) + 1.0)) / 2.0) * s;}
+ return ((y * (y * (-3.0 * y + 4.0) + 1.0)) / 2.0) * s;} //7Flops
inline float delta(float y, float s){
- return ((y * y * (y - 1.0)) / 2.0) * s;}
+ return ((y * y * (y - 1.0)) / 2.0) * s;} //5Flops
+//Total: 25Flops
diff --git a/HySoP/hysop/gpu/cl_src/weights_m6prime.cl b/HySoP/hysop/gpu/cl_src/weights_m6prime.cl
index cd21a32e74bff85d463b68da55d44bd8f83749ee..301b115d5ab1d1c8fd4ac11ffe31c739798d40f3 100644
--- a/HySoP/hysop/gpu/cl_src/weights_m6prime.cl
+++ b/HySoP/hysop/gpu/cl_src/weights_m6prime.cl
@@ -7,14 +7,15 @@
* @return weights
*/
inline float alpha(float y, float s){
- return (y * (y * (y * (y * (13.0 - 5.0 * y) - 9.0) - 1.0) + 2.0) / 24.0) * s;}
+ return (y * (y * (y * (y * (13.0 - 5.0 * y) - 9.0) - 1.0) + 2.0) / 24.0) * s;} //11Flops
inline float beta(float y, float s){
- return (y * (y * (y * (y * (25.0 * y - 64.0) + 39.0) + 16.0) - 16.0) / 24.0)*s;}
+ return (y * (y * (y * (y * (25.0 * y - 64.0) + 39.0) + 16.0) - 16.0) / 24.0)*s;} //11Flops
inline float gamma(float y, float s){
- return (y * y * (y * (y * (126.0 - 50.0 * y) - 70.0) - 30.0) / 24.0 + 1.0)*s;}
+ return (y * y * (y * (y * (126.0 - 50.0 * y) - 70.0) - 30.0) / 24.0 + 1.0)*s;} //11Flops
inline float delta(float y, float s){
- return (y * (y * (y * (y * (50.0 * y - 124.0) + 66.0) + 16.0) + 16.0) / 24.0)*s;}
+ return (y * (y * (y * (y * (50.0 * y - 124.0) + 66.0) + 16.0) + 16.0) / 24.0)*s;} //11FLops
inline float eta(float y, float s){
- return (y * (y * (y * (y * (61.0 - 25.0 * y) - 33.0) - 1.0) - 2.0) / 24.0)*s;}
+ return (y * (y * (y * (y * (61.0 - 25.0 * y) - 33.0) - 1.0) - 2.0) / 24.0)*s;} //11Flops
inline float zeta(float y, float s){
- return (y * y * y * (y * (5.0 * y - 12.0) + 7.0) / 24.0)*s;}
+ return (y * y * y * (y * (5.0 * y - 12.0) + 7.0) / 24.0)*s;} //9Flops
+//Total: 65Flops
diff --git a/HySoP/hysop/gpu/rendering.py b/HySoP/hysop/gpu/rendering.py
new file mode 100644
index 0000000000000000000000000000000000000000..aad493c9f06c839ff4862a760fe84b677639a503
--- /dev/null
+++ b/HySoP/hysop/gpu/rendering.py
@@ -0,0 +1,218 @@
+"""
+@package parmepy.gpu.rendering
+
+Handle OpenGL realtime rendering through OpenCL/OpenGL interoperability.
+@see http://enja.org/2011/03/22/adventures-in-pyopencl-part-2-particles-with-pyopengl/
+@see http://www.dyn-lab.com/articles/cl-gl.html
+"""
+from parmepy.constants import np, debug
+from parmepy.operator.monitors.monitoring import Monitoring
+from parmepy.gpu.tools import get_opengl_shared_environment, build_src
+from parmepy.fields.gpu_discrete import GPUScalarField
+from parmepy.gpu import PARMES_REAL_GPU, cl, GPU_SRC
+from parmepy.numerics.gpu_kernel import KernelLauncher
+from OpenGL.GL import *
+from OpenGL.GLU import *
+from OpenGL.GLUT import *
+
+
+class OpenGLRendering(Monitoring):
+ """
+ """
+
+ @debug
+ def __init__(self, field, delay=10, frequency=1):
+ """
+ """
+ Monitoring.__init__(self, [field], frequency)
+ print "Rendering Init"
+ if not field.dimension == 2:
+ raise ValueError("Rendering implemented in 2D only.")
+ self.delay = delay
+ #mouse handling for transforming scene
+ self.mouse_down = False
+ self.mouse_old = np.array([0., 0.])
+ self.rotate = np.array([0., 0., 0.])
+ self.translate = np.array([0., 0., 0.])
+ self.translate[0:2] = -(field.domain.origin + field.domain.length / 2.)
+ self.initrans = np.array([0., 0., -2.])
+ glutInit(sys.argv)
+ width, height = 800, 800
+ glutInitWindowSize(width, height)
+ glutInitWindowPosition(0, 0)
+ glutCreateWindow('OpenCL/OpenGL Parmepy')
+
+ glViewport(0, 0, width, height)
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(30.0, width / float(height), 1.5, 2.)
+ glMatrixMode(GL_MODELVIEW)
+
+ #handle user input
+ glutKeyboardFunc(self.on_key)
+ ## Fix mouse motion
+ #glutMouseFunc(self.on_click)
+ #glutMotionFunc(self.on_mouse_motion)
+
+ self.platform, self.device, self.ctx, self.queue = \
+ get_opengl_shared_environment(platform_id=0, device_id=0,
+ device_type='gpu',
+ precision=PARMES_REAL_GPU)
+
+ glClearColor(1, 1, 1, 1)
+ glColor(0, 0, 1)
+
+ def setUp(self):
+ print "rendering SetUp"
+ for df in self.variables[0].discreteField:
+ if isinstance(df, GPUScalarField):
+ self.gpu_field = df
+ print "Rendered field: ", self.gpu_field
+ ## Create OpenGL vbo
+ self.pos_vbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, self.pos_vbo)
+ glBufferData(GL_ARRAY_BUFFER, self.gpu_field.data.nbytes * 4,
+ None, GL_DYNAMIC_DRAW)
+ glEnableClientState(GL_VERTEX_ARRAY)
+ glVertexPointer(4, GL_FLOAT, 0, None)
+
+ self.color_vbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, self.color_vbo)
+ glBufferData(GL_ARRAY_BUFFER, self.gpu_field.data.nbytes * 4,
+ None, GL_DYNAMIC_DRAW)
+ glEnableClientState(GL_COLOR_ARRAY)
+ glColorPointer(4, GL_FLOAT, 0, None)
+
+ self.pos = cl.GLBuffer(self.ctx,
+ cl.mem_flags.READ_WRITE,
+ int(self.pos_vbo))
+ self.color = cl.GLBuffer(self.ctx,
+ cl.mem_flags.READ_WRITE,
+ int(self.color_vbo))
+ self.gl_objects = [self.pos, self.color]
+ print "Device memory used for rendering: ",
+ print self.pos.size + self.color.size, " Bytes"
+
+ f = open(GPU_SRC + 'rendering.cl')
+ gpu_src = "".join(f.readlines())
+ f.close()
+ self.prg = build_src(self.device, self.ctx, gpu_src,
+ " -cl-single-precision-constant",
+ PARMES_REAL_GPU)
+ self.numMethod = KernelLauncher(self.prg.colorize,
+ self.queue,
+ self.gpu_field.data.shape,
+ None)
+
+###GL CALLBACKS
+ def timer(self, t):
+ glutTimerFunc(t, self.timer, t)
+ glutPostRedisplay()
+
+ def on_key(self, *args):
+ ESCAPE = '\033'
+ if args[0] == ESCAPE or args[0] == 'q':
+ sys.exit()
+
+ def on_click(self, button, state, x, y):
+ if state == GLUT_DOWN:
+ self.mouse_down = True
+ self.button = button
+ else:
+ self.mouse_down = False
+ self.mouse_old[0] = x
+ self.mouse_old[1] = y
+
+ def on_mouse_motion(self, x, y):
+ dx = x - self.mouse_old[0]
+ dy = y - self.mouse_old[0]
+ if self.mouse_down and self.button == 0: # left button
+ self.rotate[0] += dy * .02
+ self.rotate[1] += dx * .02
+ elif self.mouse_down and self.button == 2: # right button
+ self.translate[2] -= dy * .001
+ self.mouse_old[0] = x
+ self.mouse_old[1] = y
+ ###END GL CALLBACKS
+
+ def setMainLoop(self, problem):
+ """
+ """
+ print "Setting rendering main loop"
+ #this will call draw every 10ms
+ glutTimerFunc(self.delay, self.timer, self.delay)
+
+ def problem_step():
+ if not problem.timer.end:
+ if (problem.domain.topologies[0].rank == 0):
+ print "==== Iteration : {0:3d} t={1:6.3f} ====".format(
+ problem.timer.ite, problem.timer.t + problem.timer.dt)
+ for op in problem.operators:
+ op.apply(problem.timer.t, problem.timer.dt,
+ problem.timer.ite)
+ problem.timer.step()
+ glutDisplayFunc(problem_step)
+
+ def finalize(self):
+ print "rendering finalization"
+
+ def apply(self, t, dt, ite):
+ """
+ """
+ self.queue.finish()
+ self.numMethod(self.gpu_field.gpu_data, self.pos, self.color)
+ cl.enqueue_release_gl_objects(self.queue, self.gl_objects)
+ self.queue.finish()
+
+ glFlush()
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
+ glMatrixMode(GL_MODELVIEW)
+ glLoadIdentity()
+
+ #handle mouse transformations
+ glTranslatef(self.initrans[0], self.initrans[1], self.initrans[2])
+ glRotatef(self.rotate[0], 1, 0, 0)
+ glRotatef(self.rotate[1], 0, 1, 0) # we switched around the axis so make this rotate_z
+ glTranslatef(self.translate[0], self.translate[1], self.translate[2])
+
+ glEnable(GL_POINT_SMOOTH)
+ glPointSize(3)
+ #glEnable(GL_BLEND)
+ glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
+
+ glBindBuffer(GL_ARRAY_BUFFER, self.color_vbo)
+ glColorPointer(4, GL_FLOAT, 0, None)
+ glBindBuffer(GL_ARRAY_BUFFER, self.pos_vbo)
+ glVertexPointer(4, GL_FLOAT, 0, None)
+
+ glEnableClientState(GL_VERTEX_ARRAY)
+ glEnableClientState(GL_COLOR_ARRAY)
+ glDrawArrays(GL_POINTS, 0, np.prod(self.gpu_field.data.shape))
+ glDisableClientState(GL_COLOR_ARRAY)
+ glDisableClientState(GL_VERTEX_ARRAY)
+
+ #glDisable(GL_BLEND)
+ orig = np.array([0, 0, 0])
+ x_one = np.array([1, 0, 0])
+ y_one = np.array([0, 1, 0])
+ z_one = np.array([0, 0, 1])
+ glColor3f(1, 0, 0) # red
+ self.draw_line(orig, x_one)
+ glColor3f(0, 1, 0) # green
+ self.draw_line(orig, y_one)
+ #glColor3f(0, 0, 1) # blue
+ #self.draw_line(orig, z_one)
+
+ glFlush()
+ cl.enqueue_acquire_gl_objects(self.queue, self.gl_objects)
+
+ def draw_line(self, v1, v2):
+ glBegin(GL_LINES)
+ glVertex3f(v1[0], v1[1], v1[2])
+ glVertex3f(v2[0], v2[1], v2[2])
+ glEnd()
+
+ def startMainLoop(self):
+ """
+ """
+ glutMainLoop()
diff --git a/HySoP/hysop/gpu/tests/test_opencl_environment.py b/HySoP/hysop/gpu/tests/test_opencl_environment.py
new file mode 100644
index 0000000000000000000000000000000000000000..1e107c4eed46df1966922e6306546c9e523c3bfa
--- /dev/null
+++ b/HySoP/hysop/gpu/tests/test_opencl_environment.py
@@ -0,0 +1,52 @@
+from parmepy.gpu.tools import get_opencl_environment, \
+ get_opengl_shared_environment
+from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU
+
+
+def test_queue_unique_creation():
+ """
+ Testing that only one queue is created when multiples calls to get
+ an environment.
+ """
+ platform, device, ctx, queue = get_opencl_environment(0, 0, 'gpu',
+ PARMES_REAL_GPU)
+ platform_id = id(platform)
+ device_id = id(device)
+ ctx_id = id(ctx)
+ queue_id = id(queue)
+ platformb, deviceb, ctxb, queueb = get_opencl_environment(0, 0, 'gpu',
+ PARMES_REAL_GPU)
+ platformb_id = id(platformb)
+ deviceb_id = id(deviceb)
+ ctxb_id = id(ctxb)
+ queueb_id = id(queueb)
+ assert platformb_id == platform_id
+ assert deviceb_id == device_id
+ assert ctxb_id == ctx_id
+ assert queueb_id == queue_id
+
+
+def test_queue_with_gl_sharing():
+ """
+ Testing that queue and context is changed when OpenGL rendering is
+ activated. Environment, platform and device remains unchanged.
+ """
+ from OpenGL.GLUT import glutInit
+ import sys
+ glutInit(sys.argv)
+ platform, device, ctx, queue = get_opengl_shared_environment(
+ 0, 0, 'gpu', PARMES_REAL_GPU)
+ platform_id = id(platform)
+ device_id = id(device)
+ ctx_id = id(ctx)
+ queue_id = id(queue)
+ platformb, deviceb, ctxb, queueb = get_opencl_environment(0, 0, 'gpu',
+ PARMES_REAL_GPU)
+ platformb_id = id(platformb)
+ deviceb_id = id(deviceb)
+ ctxb_id = id(ctxb)
+ queueb_id = id(queueb)
+ assert platformb_id == platform_id
+ assert deviceb_id == device_id
+ assert ctxb_id == ctx_id
+ assert queueb_id == queue_id
diff --git a/HySoP/hysop/gpu/tools.py b/HySoP/hysop/gpu/tools.py
index 7cd2711a88c787f62207e7dc36eacf8bd36dbbb3..4c39c1281489ca632a4c3171489595bf5a83ef32 100644
--- a/HySoP/hysop/gpu/tools.py
+++ b/HySoP/hysop/gpu/tools.py
@@ -3,66 +3,145 @@
Tools for gpu management.
"""
-from parmepy import __VERBOSE__
+#from parmepy import __VERBOSE__
from parmepy.gpu import cl, PARMES_REAL_GPU
+__VERBOSE__ = True
__platform, __device, __ctx, __queue = None, None, None, None
__precision = None
+def _get_platform(platform_id):
+ """
+ """
+ #Get platform.
+ try:
+ ## OpenCL platform
+ platform = cl.get_platforms()[platform_id]
+ except IndexError:
+ print " Incorrect platform_id :", platform_id, ".",
+ print " Only ", len(cl.get_platforms()), " available.",
+ print " Getting defalut platform. "
+ __platform = cl.get_platforms()[0]
+ if __VERBOSE__:
+ print " Platform "
+ print " - Name :", platform.name
+ print " - Version :", platform.version
+ return platform
+
+
+def _get_device(platform, device_id, device_type):
+ """
+ """
+ try:
+ ## OpenCL device
+ device = platform.get_devices(
+ eval("cl.device_type." + device_type.upper()))[device_id]
+ except cl.RuntimeError as e:
+ print "RuntimeError:", e
+ device = cl.create_some_context().devices[0]
+ except AttributeError as e:
+ print "AttributeError:", e
+ device = cl.create_some_context().devices[0]
+ if __VERBOSE__:
+ print " Device"
+ print " - Name :",
+ print device.name
+ print " - Type :",
+ print cl.device_type.to_string(device.type)
+ print " - C Version :",
+ print device.opencl_c_version
+ print " - Global mem size :",
+ print device.global_mem_size / (1024 ** 3), "GB"
+ return device
+
+
+def _get_context(device, props=None):
+ """
+ """
+ if props is None:
+ ctx = cl.Context([device])
+ else:
+ ctx = cl.Context(properties=props, devices=[device])
+ if __VERBOSE__:
+ print " Context:"
+ print " - properties :", props
+ return ctx
+
+
+def _get_queue(ctx):
+ """
+ """
+ props = cl.command_queue_properties.PROFILING_ENABLE
+ queue = cl.CommandQueue(ctx, properties=props)
+ if __VERBOSE__:
+ print " Queue"
+ print " - properties :", props
+ print "==="
+ return queue
+
+
+def get_opengl_shared_environment(platform_id, device_id, device_type,
+ precision):
+ """
+ """
+ global __platform
+ global __device
+ global __ctx
+ global __queue
+ global __precision
+ print "Create OpenGL,", __platform, __device, __ctx, __queue
+ if __precision is None and not precision is None:
+ __precision = precision
+ if __queue is None:
+ if __VERBOSE__:
+ print "=== OpenGL shared environment ==="
+ __platform = _get_platform(platform_id)
+ #Get device.
+ __device = _get_device(__platform, device_id, device_type)
+ #Creates GPU Context
+ ## OpenCL context
+ from pyopencl.tools import get_gl_sharing_context_properties
+ import sys
+ if sys.platform == "darwin":
+ properties = get_gl_sharing_context_properties()
+ else:
+ # Some OSs prefer clCreateContextFromType, some prefer
+ # clCreateContext. Try both.
+ try:
+ properties = [
+ (cl.context_properties.PLATFORM, platform)]
+ except:
+ properties = [
+ (cl.context_properties.PLATFORM, platform)]
+ __ctx = _get_context(__device, properties)
+ #Create CommandQueue on the GPU Context
+ ## OpenCL command queue
+ __queue = _get_queue(__ctx)
+ return __platform, __device, __ctx, __queue
+
+
def get_opencl_environment(platform_id, device_id, device_type, precision):
global __platform
global __device
global __ctx
global __queue
global __precision
+ print "Create OpenCL,", __platform, __device, __ctx, __queue
if __precision is None and not precision is None:
__precision = precision
if __queue is None:
if __VERBOSE__:
print "=== OpenCL environment ==="
- #Get platform.
- try:
- ## OpenCL platform
- __platform = cl.get_platforms()[platform_id]
- except IndexError:
- print " Incorrect platform_id :", platform_id, ".",
- print " Only ", len(cl.get_platforms()), " available.",
- print " Getting defalut platform. "
- __platform = cl.get_platforms()[0]
- if __VERBOSE__:
- print " Platform "
- print " - Name :", __platform.name
- print " - Version :", __platform.version
+ __platform = _get_platform(platform_id)
#Get device.
- try:
- ## OpenCL device
- __device = __platform.get_devices(
- eval("cl.device_type." + device_type.upper()))[device_id]
- except cl.RuntimeError as e:
- print "RuntimeError:", e
- __device = cl.create_some_context().devices[0]
- except AttributeError as e:
- print "AttributeError:", e
- __device = cl.create_some_context().devices[0]
- if __VERBOSE__:
- print " Device"
- print " - Name :",
- print __device.name
- print " - Type :",
- print cl.device_type.to_string(__device.type)
- print " - C Version :",
- print __device.opencl_c_version
- print " - Global mem size :",
- print __device.global_mem_size / (1024 ** 3), "GB"
- print "===\n"
+ __device = _get_device(__platform, device_id, device_type)
#Creates GPU Context
## OpenCL context
- __ctx = cl.Context([__device])
+ __ctx = _get_context(__device)
#Create CommandQueue on the GPU Context
## OpenCL command queue
- __queue = cl.CommandQueue(
- __ctx, properties=cl.command_queue_properties.PROFILING_ENABLE)
+ __queue = _get_queue(__ctx)
return __platform, __device, __ctx, __queue
@@ -138,13 +217,32 @@ def explore():
print "OpenCL exploration : "
platforms = cl.get_platforms()
platforms_info = ["name", "version", "vendor", "profile", "extensions"]
- devices_info = ["name", "version", "vendor", "profile", "extensions",
- "available", "type", "compiler_available", "double_fp_config", "single_fp_config",
- "global_mem_size", "global_mem_cache_type", "global_mem_cache_size", "global_mem_cacheline_size",
- "local_mem_size", "local_mem_type", "max_clock_frequency", "max_compute_units",
- "max_constant_buffer_size", "max_mem_alloc_size", "max_work_group_size",
- "max_work_item_dimensions", "max_work_item_sizes",
- "preferred_vector_width_double", "preferred_vector_width_float", "preferred_vector_width_int"]
+ devices_info = ["name",
+ "version",
+ "vendor",
+ "profile",
+ "extensions",
+ "available",
+ "type",
+ "compiler_available",
+ "double_fp_config",
+ "single_fp_config",
+ "global_mem_size",
+ "global_mem_cache_type",
+ "global_mem_cache_size",
+ "global_mem_cacheline_size",
+ "local_mem_size",
+ "local_mem_type",
+ "max_clock_frequency",
+ "max_compute_units",
+ "max_constant_buffer_size",
+ "max_mem_alloc_size",
+ "max_work_group_size",
+ "max_work_item_dimensions",
+ "max_work_item_sizes",
+ "preferred_vector_width_double",
+ "preferred_vector_width_float",
+ "preferred_vector_width_int"]
for pltfm in platforms:
print "Platform:", pltfm.name
for pltfm_info in platforms_info:
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index b0b5f87301807e167e5e60cadaab4797b8db3a01..97de7046090094f3092bced4b8e401543f7562ce 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -106,7 +106,7 @@ class Advection(Operator):
**self.config)
else:
print "Using default advection operator"
- if self.advectedFields.vector:
+ if self.advectedFields.isVector:
raise ValueError("Not implemented yet. default" +
" advection is only for " +
"scalar fields")
diff --git a/HySoP/hysop/operator/analytic.py b/HySoP/hysop/operator/analytic.py
index a4c2f3733472516e14e1baece11c9e76de7864f7..7074f37961acbe4b4987f4896e552f9e31200109 100644
--- a/HySoP/hysop/operator/analytic.py
+++ b/HySoP/hysop/operator/analytic.py
@@ -25,7 +25,7 @@ class Analytic(Operator):
@param fields : field
"""
- if not callable(formula):
+ if not callable(formula) and 'gpu' in other_config.keys():
raise ValueError("Analytic formula must be a function")
Operator.__init__(self, [fields])
## velocity variable
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 43d7b084548ced0c929a237e412f31740a462d78..83a2cdc6c94c6012c31ed1e92e7c4ca65f8f94da 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -47,7 +47,7 @@ class Printer(Monitoring):
res = {}
for f in self.variables:
for df in f.discreteField:
- if f.vector:
+ if f.isVector:
for d in xrange(len(df.data)):
if len(df.data[d].shape) == 2:
res[df.name + S_DIR[d]] = np.expand_dims(
@@ -106,7 +106,7 @@ class Printer(Monitoring):
f.write("{0:8.12} {1:8.12} ".format(
coords[0][i, 0], coords[1][0, j]))
for field in self.variables:
- if field.vector:
+ if field.isVector:
f.write("{0:8.12} {1:8.12} ".format(
field.discreteField[0][0][i, j],
field.discreteField[0][1][i, j]))
@@ -123,7 +123,7 @@ class Printer(Monitoring):
coords[1][0, j, 0],
coords[2][0, 0, k]))
for field in self.variables:
- if field.vector:
+ if field.isVector:
f.write("{0:8.12} {1:8.12} {2:8.12} ".format(
field.discreteField[0][0][i, j, k],
field.discreteField[0][1][i, j, k],
diff --git a/HySoP/hysop/operator/particle_advection.py b/HySoP/hysop/operator/particle_advection.py
index 7a2366f95ad1004b2e4e4471b52a14fbf6c7bbc7..bd933f39d7ab1a1f12b3e2f1b6b63ab73daa34eb 100644
--- a/HySoP/hysop/operator/particle_advection.py
+++ b/HySoP/hysop/operator/particle_advection.py
@@ -56,16 +56,14 @@ class ParticleAdvection(DiscreteOperator):
## Result position
particle_position = Field(self.advec.advectedFields.domain,
"Particle_Position",
- vector=False)
- particle_position.discretize(
- self.advec.resolutions[self.advec.advectedFields])
+ isVector=False)
+ particle_position.discretize(self.advectedFields.topology)
self.part_position = particle_position.discreteField[0]
## Result advectedFields
particle_advectedFields = Field(self.advec.advectedFields.domain,
"Particle_AdvectedFields",
- vector=False)
- particle_advectedFields.discretize(
- self.advec.resolutions[self.advec.advectedFields])
+ isVector=False)
+ particle_advectedFields.discretize(self.advectedFields.topology)
self.part_advectedFields = particle_advectedFields.discreteField[0]
if self.method.find('rk2') >= 0:
self.name += '_rk2'
diff --git a/HySoP/hysop/operator/tests/test_particle_advection.py b/HySoP/hysop/operator/tests/test_particle_advection.py
index 6706b7599c4aa6f0482d6d41f0398d7cd90ec50a..ec8bba350018215c2c18e117bc115108120ccbf7 100644
--- a/HySoP/hysop/operator/tests/test_particle_advection.py
+++ b/HySoP/hysop/operator/tests/test_particle_advection.py
@@ -11,11 +11,12 @@ def test_null_velocity():
dom = Box(dimension=2, origin=[0., 0.], length=[1., 1.])
resolTopo = [33, 33]
csf = AnalyticalField(dom, formula=lambda x, y: 1.)
- cvf = AnalyticalField(dom, vector=True, formula=lambda x, y: (0., 0.))
- advec = Advection(cvf, csf)
- advec.setUp(resolutions={csf: [33, 33],
- cvf: [33, 33]},
- method='')
+ cvf = AnalyticalField(dom, isVector=True, formula=lambda x, y: (0., 0.))
+ advec = Advection(cvf, csf,
+ resolutions={csf: resolTopo,
+ cvf: resolTopo},
+ method='')
+ advec.setUp()
pytest.raises(ValueError, advec.apply, 0., 1., 1)
# advec.apply(0., 1., 1)
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 961d298cc796e1c97c34ded2a54e937fe6180540..5fcb4bcbdb1434131399f738a3cf06f31fc4d06c 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -7,18 +7,26 @@ from abc import ABCMeta, abstractmethod
from parmepy.constants import debug
from parmepy import __VERBOSE__
from parmepy.operator.monitors.monitoring import Monitoring
+try:
+ from parmepy.gpu.rendering import OpenGLRendering
+except ImportError:
+ OpenGLRendering = None
class Problem(object):
"""
Problem representation.
- Contains several operators that apply on variables.
+ Contains several operators and monitors that apply on variables.
Variables are defined on different domains.\n
- The problem contains a solver that define how the problem is solved.\n
- A problem is based on a topology that deploys
- its resolution over MPI processus.\n
- Finally, the problem handle outputs thanks to a io object.
+ Each operator is set up and variables are initialized in a set up step.\n
+ To solve the problem, a loop over time-steps is launched. A step consists
+ in calling the apply method of each operators and monitors.\n
+ To finish, a finalize method is called.\
+
+ \remark: In case of GPU real-time rendering (i.e. use of an
+ OpenGLRendering object), The loop over time-steps is passed to the
+ OpenGL Utility Toolkit (GLUT).
"""
__metaclass__ = ABCMeta
@@ -34,6 +42,7 @@ class Problem(object):
Create a transport problem instance.
@param operators : list of operators.
+ @param monitors : list of monitors.
"""
## Problem operators
self.operators = operators
@@ -66,6 +75,13 @@ class Problem(object):
v.initialize()
if __VERBOSE__:
print "===="
+ ## Set the callback function to OpenGL
+ self.glRenderer = None
+ if not OpenGLRendering is None:
+ for op in self.operators:
+ if isinstance(op, OpenGLRendering):
+ self.glRenderer = op
+ op.setMainLoop(self)
@debug
def solve(self):
@@ -77,20 +93,22 @@ class Problem(object):
At end of time step, call an io step.\n
Displays timings at simulation end.
"""
- ## Iterations counter
- ite = 0
print "\n\n Start solving ..."
for op in self.operators:
if isinstance(op, Monitoring):
- op.apply(self.timer.t, self.timer.dt, ite)
- while not self.timer.end:
- ite = ite + 1
- if (self.domain.topologies[0].rank == 0):
- print "==== Iteration : {0:3d} t={1:6.3f} ====".format(
- self.timer.ite, self.timer.t + self.timer.dt)
- for op in self.operators:
- op.apply(self.timer.t, self.timer.dt, ite)
- self.timer.step()
+ op.apply(self.timer.t, self.timer.dt, ite=0)
+ ## Pass the main loop to GLUT in case of GL_Rendering
+ if not self.glRenderer is None:
+ self.glRenderer.startMainLoop()
+ ## Or run a loop normally
+ else:
+ while not self.timer.end:
+ if (self.domain.topologies[0].rank == 0):
+ print "==== Iteration : {0:3d} t={1:6.3f} ====".format(
+ self.timer.ite, self.timer.t + self.timer.dt)
+ for op in self.operators:
+ op.apply(self.timer.t, self.timer.dt, self.timer.ite)
+ self.timer.step()
print "\n\n End solving\n"
print "==== Finalization ===="
for op in self.operators:
@@ -146,7 +164,7 @@ class Timer:
## Is simulation is terminated
self.end = False
## Iteration counter
- self.ite = 0
+ self.ite = 1
def step(self):
"""