From 79c4a326d636c53d1e97a3b7b8d48b53b5ff93bb Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Wed, 11 Jul 2012 18:01:07 +0000
Subject: [PATCH] Adding GPU solver (initialization only)
---
Examples/mainJM.py | 15 ++++---
Examples/mainJM_kernels.cl | 37 +++++++++--------
HySoP/hysop/__init__.py.in | 28 +++++++++++--
HySoP/hysop/constants.py | 2 +
HySoP/hysop/fields/analytical.py | 2 +
HySoP/hysop/fields/continuous.py | 10 +++--
HySoP/hysop/fields/discrete.py | 41 +++++++++++++------
.../particular_solvers/GPUParticularSolver.py | 2 +-
.../GPUParticularSolver_GLRender.py | 2 +-
HySoP/hysop/particular_solvers/__init__.py | 6 +++
HySoP/hysop/particular_solvers/basic.py | 5 ++-
HySoP/hysop/problem/problem.py | 11 ++++-
HySoP/hysop/tools/__init__.py | 6 +++
HySoP/setup.py.in | 11 +++--
HySoP/src/main/main.f90 | 2 -
15 files changed, 125 insertions(+), 55 deletions(-)
diff --git a/Examples/mainJM.py b/Examples/mainJM.py
index 4723fbfec..8aa4df466 100644
--- a/Examples/mainJM.py
+++ b/Examples/mainJM.py
@@ -20,7 +20,7 @@ def scalaire(x, y, z):
def run():
# Parameters
- timeStep = 0.02
+ timeStep = 1. # 0.02
finalTime = 1.
outputFilePrefix = './res/RK2_'
outputModulo = 0
@@ -31,20 +31,23 @@ def run():
box = pp.Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
## Fields
- scal = pp.AnalyticalField(domain=box, formula=scalaire, name='Scalar')
- velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
+ #scal = pp.AnalyticalField(domain=box, formula=scalaire, name='Scalar')
+ #velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
+ scal = pp.ContinuousField(domain=box, name='Scalar')
+ velo = pp.ContinuousField(domain=box, name='Velocity', vector=True)
## Operators
advec = pp.Transport(velo, scal)
## Solver creation (discretisation of objects is done in solver initialisation)
- topo3D = pp.CartesianTopology(domain=box, resolution=[4, 4, 4], dim=3)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[256, 256, 256], dim=3)
##Problem
pb = pp.Problem(topo3D, [advec])
## Setting solver to Problem
- pb.setSolver(pp.ParticularSolver, finalTime, timeStep)
+ pb.setSolver(finalTime, timeStep, 'gpu', src='./mainJM_kernels.cl')
+ #pb.setSolver(finalTime, timeStep, 'gpu')
# ## Setting io to Problem
pb.setIO(pp.Printer(fields=[scal, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
@@ -52,7 +55,7 @@ def run():
t1 = time.time()
## Solve problem
- pb.solve()
+ #pb.solve()
tf = time.time()
diff --git a/Examples/mainJM_kernels.cl b/Examples/mainJM_kernels.cl
index 007e7c297..a45672ba4 100644
--- a/Examples/mainJM_kernels.cl
+++ b/Examples/mainJM_kernels.cl
@@ -1,16 +1,17 @@
__kernel void initScalar(__global float* values,
- __private const float t,
__private const float4 min,
- __private const uint4 nb,
__private const float4 size)
{
- __private uint ind,ix,iy,iz;
+ __private uint ind,ix,iy,iz,nbx,nby,nbz;
__private float px,py,pz,s;
ix = get_global_id(0);
iy = get_global_id(1);
iz = get_global_id(2);
+ nbx = get_global_size(0);
+ nby = get_global_size(1);
+ nbz = get_global_size(2);
px = min.x + (float)(ix)*size.x;
py = min.y + (float)(iy)*size.y;
@@ -21,32 +22,34 @@ __kernel void initScalar(__global float* values,
else
s = 0.0f;
// Write
- ind = iy + ix*nb.y;
+ ind = iz*nby*nbx + iy*nbx + ix;
values[ind] = s;
}
// velocity field
-__kernel void initVelocity(__global float* values,
- __private const float t,
- __private const float4 min,
- __private const uint4 nb,
- __private const float4 size
- )
+__kernel void initVelocity(__global float* values_x,
+ __global float* values_y,
+ __global float* values_z,
+ __private const float4 min,
+ __private const float4 size
+ )
{
- __private uint ix,iy,iz, ind;
+ __private uint ix,iy,iz,nbx,nby,nbz, ind;
ix = get_global_id(0);
iy = get_global_id(1);
iz = get_global_id(2);
+ nbx = get_global_size(0);
+ nby = get_global_size(1);
+ nbz = get_global_size(2);
+
// Write x component
- ind = ix*nb.y*nb.z + iy*nb.z + iz;
- values[ind] = 1.0f;
+ ind = iz*nby*nbx + iy*nbx + ix;
+ values_x[ind] = 1.0f;
// Write y component
- ind = (nb.x*nb.y*nb.z) + iy*nb.x*nb.z + ix*nb.z + iz;
- values[ind] = 1.0f;
+ values_y[ind] = 1.0f;
// Write z component
- ind = (nb.x*nb.y*nb.z)*2 + iz*nb.x*nb.y + ix*nb.y + iy;
- values[ind] = 1.0f;
+ values_z[ind] = 1.0f;
}
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index 4e4aad24f..85cbbdae5 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -4,8 +4,8 @@
Python package dedicated to flow simulation using particular methods on hybrid architectures (MPI-GPU)
"""
-__version__=1.00
-__all__=['Box','CartesianTopology','ScalarField']
+__version__ = 1.00
+__all__ = ['Box', 'CartesianTopology', 'ScalarField']
import os
import site
@@ -13,12 +13,12 @@ import mpi4py.MPI as MPI
rank_world = MPI.COMM_WORLD.Get_rank()
if(rank_world == 0):
- print "Starting @PACKAGE_NAME@ version "+str(__version__)+".\n"
+ print "Starting @PACKAGE_NAME@ version " + str(__version__) + ".\n"
import constants
import domain.box
-## Box-type physical domain
+## Box-type physical domain
Box = domain.box.Box
## Cartesian grid
@@ -30,5 +30,25 @@ CartesianTopology = fields.topology.CartesianTopology
## Fields
import fields.discrete
import fields.continuous
+import fields.analytical
ScalarField = fields.discrete.ScalarField
ContinuousField = fields.continuous.ContinuousField
+AnalyticalField = fields.analytical.AnalyticalField
+
+## Operators
+import operator.Transport
+Transport = operator.Transport.Transport
+
+## Problem
+import problem.problem
+Problem = problem.problem.Problem
+
+## Solver
+import ParticularSolvers.basic
+import ParticularSolvers.gpu
+ParticleSolver = ParticularSolvers.basic.ParticleSolver
+GPUParticleSolver = ParticularSolvers.gpu.GPUParticleSolver
+
+## Tools
+import tools.printer
+Printer = tools.printer.Printer
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
index 1af66a9bc..6b4b94b95 100644
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -27,3 +27,5 @@ ZDIR = 2
PERIODIC = 0
## Directions string
S_DIR = ["_X", "_Y", "_Z"]
+## GPU deflault sources
+GPU_SRC = "/Users/jmetancelin/These/SVN/parmes/Parmes/parmepy/ParticularSolvers/gpu_src.cl"
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
index d4ec1738f..2b1dcc353 100644
--- a/HySoP/hysop/fields/analytical.py
+++ b/HySoP/hysop/fields/analytical.py
@@ -39,8 +39,10 @@ class AnalyticalField(ContinuousField):
def initialize(self):
if self._fieldId == -1:
raise ValueError("Cannot initialise analytical field non discretized.")
+ print self.name,
for dF in self.discreteField:
dF.initialize(self.formula)
+ print "Done"
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 7bfb54eb6..aa088c6a2 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -39,21 +39,23 @@ class ContinuousField(object):
"""
self._fieldId += 1
- print "Discretisation de ", self.name,
if self.vector:
- print "Vector Field"
self.discreteField.append(VectorField(self, topology,
name=self.name + "_D_" + str(self._fieldId),
idFromParent=self._fieldId))
else:
- print "Scalar Field"
self.discreteField.append(ScalarField(self, topology,
name=self.name + "_D_" + str(self._fieldId),
idFromParent=self._fieldId))
return self.discreteField[self._fieldId], self._fieldId
def initialize(self):
- pass
+ if self._fieldId == -1:
+ raise ValueError("Cannot initialise analytical field non discretized.")
+ print self.name,
+ for dF in self.discreteField:
+ dF.initialize()
+ print "Done"
def __str__(self):
""" Field info display """
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
index bccf1be47..f5f77293d 100644
--- a/HySoP/hysop/fields/discrete.py
+++ b/HySoP/hysop/fields/discrete.py
@@ -22,12 +22,15 @@ class ScalarField(object):
self.name = name
self.dimension = topology.domain.dimension
- resolution = self.topology.mesh.resolution
+ self.resolution = self.topology.mesh.resolution
## \todo is numpy array zeros or empty???
- self.data = np.zeros((resolution), dtype=PARMES_REAL, order=ORDER)
+ self.data = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
+ self.gpu_data = None
self.__parentField = parent
self.__idFromParent = idFromParent
self.domain = self.__parentField.domain.discreteDomain
+ self.contains_data = False
+ self.vector = False
def __str__(self):
""" Display class information """
@@ -54,16 +57,21 @@ class ScalarField(object):
"""
Initialize values with given formula
"""
- v_formula = np.vectorize(formula)
- if self.dimension == 3:
- self.data = v_formula(self.topology.mesh.coords[0],
- self.topology.mesh.coords[1],
- self.topology.mesh.coords[2])
- elif self.dimension == 2:
- self.data = v_formula(self.topology.mesh.coords[0],
- self.topology.mesh.coords[1])
+ if formula is not None:
+ print "...",
+ v_formula = np.vectorize(formula)
+ if self.dimension == 3:
+ self.data = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1],
+ self.topology.mesh.coords[2])
+ elif self.dimension == 2:
+ self.data = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1])
+ else:
+ self.data = v_formula(self.topology.mesh.coords[0])
+ self.contains_data = True
else:
- self.data = v_formula(self.topology.mesh.coords[0])
+ print "",
def tofile(self, filename):
evtk.imageToVTK(filename, pointData={self.name: self.data})
@@ -83,12 +91,15 @@ class VectorField(object):
self.name = name
self.dimension = topology.domain.dimension
- resolution = self.topology.mesh.resolution
+ self.resolution = self.topology.mesh.resolution
## \todo is numpy array zeros or empty???
- self.data = [np.zeros((resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dimension)]
+ self.data = [np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dimension)]
+ self.gpu_data = [None for d in xrange(self.dimension)]
self.__parentField = parent
self.__idFromParent = idFromParent
self.domain = self.__parentField.domain.discreteDomain
+ self.contains_data = False
+ self.vector = True
def __str__(self):
""" Display class information """
@@ -137,6 +148,7 @@ class VectorField(object):
"""
#print self.topology.mesh.coords
if formula is not None:
+ print "...",
v_formula = np.vectorize(formula)
if self.dimension == 3:
self.data[0], self.data[1], self.data[2] = v_formula(self.topology.mesh.coords[0],
@@ -145,6 +157,9 @@ class VectorField(object):
elif self.dimension == 2:
self.data[0], self.data[1] = v_formula(self.topology.mesh.coords[0],
self.topology.mesh.coords[1])
+ self.contains_data = True
+ else:
+ print "No formula",
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/particular_solvers/GPUParticularSolver.py b/HySoP/hysop/particular_solvers/GPUParticularSolver.py
index a5f6cbc64..f87f51852 100644
--- a/HySoP/hysop/particular_solvers/GPUParticularSolver.py
+++ b/HySoP/hysop/particular_solvers/GPUParticularSolver.py
@@ -6,7 +6,7 @@ GPU Particular solver description.
import os
import pyopencl as cl
from ..Param import *
-from particular_solver import ParticularSolver
+from ParticularSolver import ParticularSolver
from ..Domain.ParticleField import ParticleField
from ..Variable.DiscreteVariable import DiscreteVariable
from ..Operator.RemeshingDOp import RemeshingDOp
diff --git a/HySoP/hysop/particular_solvers/GPUParticularSolver_GLRender.py b/HySoP/hysop/particular_solvers/GPUParticularSolver_GLRender.py
index 90a16bf04..24432e771 100644
--- a/HySoP/hysop/particular_solvers/GPUParticularSolver_GLRender.py
+++ b/HySoP/hysop/particular_solvers/GPUParticularSolver_GLRender.py
@@ -7,7 +7,7 @@ import os
import sys
import pyopencl as cl
from ..Param import *
-from particular_solver import ParticularSolver
+from ParticularSolver import ParticularSolver
from GPUParticularSolver import GPUParticularSolver
from ..Domain.ParticleField import ParticleField
from ..Variable.DiscreteVariable import DiscreteVariable
diff --git a/HySoP/hysop/particular_solvers/__init__.py b/HySoP/hysop/particular_solvers/__init__.py
index e69de29bb..6324d8802 100644
--- a/HySoP/hysop/particular_solvers/__init__.py
+++ b/HySoP/hysop/particular_solvers/__init__.py
@@ -0,0 +1,6 @@
+"""
+@package parmepy.ParticularSolvers
+
+Everything concerning particle solvers.
+
+"""
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index 16ffa8a14..c5783fd23 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -12,7 +12,7 @@ from ..tools.timer import Timer
from ..tools.printer import Printer
-class ParticularSolver(Solver):
+class ParticleSolver(Solver):
"""
Particular solver description.
Link with differents numericals methods used.
@@ -68,7 +68,10 @@ class ParticularSolver(Solver):
## Discretise fields and initialize
for v in self.problem.variables:
v.discretize(self.problem.topology)
+ print "=== Fields initialization ==="
+ for v in self.problem.variables:
v.initialize()
+ print "==\n"
## Discretise operators
for op in self.problem.operators:
if op is self.advection:
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 74e40d065..fb209fcd7 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -3,6 +3,8 @@
@package Problem
Problem representation
"""
+from ..particular_solvers.basic import ParticleSolver
+from ..particular_solvers.gpu import GPUParticleSolver
class Problem():
@@ -34,8 +36,13 @@ class Problem():
def setIO(self, io):
self.io = io
- def setSolver(self, solver, t_end, dt):
- self.solver = solver(self, t_end, dt)
+ def setSolver(self, t_end, dt, solver_type='basic', **kwargs):
+ if solver_type == 'basic':
+ self.solver = ParticleSolver(self, t_end, dt, **kwargs)
+ elif solver_type == 'gpu':
+ self.solver = GPUParticleSolver(self, t_end, dt, **kwargs)
+ else:
+ raise ValueError("Unknown solver type : " + str(solver_type))
self.solver.initialize()
def solve(self):
diff --git a/HySoP/hysop/tools/__init__.py b/HySoP/hysop/tools/__init__.py
index e69de29bb..44a793122 100644
--- a/HySoP/hysop/tools/__init__.py
+++ b/HySoP/hysop/tools/__init__.py
@@ -0,0 +1,6 @@
+"""
+@package parmepy.tools
+
+Everything concerning tools.
+
+"""
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index e74a6d9d0..17d235b73 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -37,8 +37,6 @@ parmes_libdir = ['@CMAKE_BINARY_DIR@/src']
parmeslib = ['@PARMES_LIBRARY_NAME@']
f2py_options = ['--no-lower']
-
-
scales_src = glob.glob('@CMAKE_SOURCE_DIR@/parmepy/scales2py/*.f90')
scalesModule = Extension(name='parmepy.scales2py',
f2py_options=f2py_options,
@@ -51,8 +49,13 @@ scalesModule = Extension(name='parmepy.scales2py',
ext_modules = [scalesModule]
test_src = glob.glob('@CMAKE_SOURCE_DIR@/parmepy/testf2py/*.f90')
-testf2pyModule=Extension(name='parmepy.testf2py',f2py_options=f2py_options,sources=test_src,include_dirs=parmes_dir,
- library_dirs=parmes_libdir,libraries=parmeslib,define_macros = [('F2PY_REPORT_ON_ARRAY_COPY','1')])
+testf2pyModule = Extension(name='parmepy.testf2py',
+ f2py_options=f2py_options,
+ sources=test_src,
+ include_dirs=parmes_dir,
+ library_dirs=parmes_libdir,
+ libraries=parmeslib,
+ define_macros=[('F2PY_REPORT_ON_ARRAY_COPY', '1')])
ext_modules.append(testf2pyModule)
diff --git a/HySoP/src/main/main.f90 b/HySoP/src/main/main.f90
index 2eca3a4d4..2d75a8941 100755
--- a/HySoP/src/main/main.f90
+++ b/HySoP/src/main/main.f90
@@ -3,7 +3,6 @@ program mainParmes
use client_data
use advec
-use solverpoisson
implicit none
@@ -17,7 +16,6 @@ call MPI_COMM_SIZE(MPI_COMM_WORLD,nbprocs,info)
print *, "hello ..."
-call testpoisson()
!!call toto(4)
call MPI_Finalize(info)
--
GitLab