From d9e4e036887a024649e5dc3072315e6cee60efa8 Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Fri, 31 Aug 2012 10:12:46 +0000
Subject: [PATCH] Add analytical Velocity time dependant
---
Examples/levelSet2D.cl | 12 ++--
Examples/levelSet2D.py | 7 ++-
Examples/levelSet3D.cl | 25 ++++----
Examples/levelSet3D.py | 5 +-
HySoP/hysop/__init__.py.in | 2 +
HySoP/hysop/operator/Velocity.py | 52 +++++++++++++++++
HySoP/hysop/operator/advection_d.py | 1 +
HySoP/hysop/operator/remeshing.py | 1 +
HySoP/hysop/operator/splitting.py | 1 +
HySoP/hysop/operator/velocity_d.py | 78 +++++++++++++++++++++++++
HySoP/hysop/particular_solvers/basic.py | 4 +-
HySoP/hysop/particular_solvers/gpu.py | 19 ++++--
12 files changed, 178 insertions(+), 29 deletions(-)
create mode 100644 HySoP/hysop/operator/Velocity.py
create mode 100644 HySoP/hysop/operator/velocity_d.py
diff --git a/Examples/levelSet2D.cl b/Examples/levelSet2D.cl
index 6f48b9e74..49e711316 100644
--- a/Examples/levelSet2D.cl
+++ b/Examples/levelSet2D.cl
@@ -13,17 +13,21 @@ __kernel void initScalar(__global float* scalar,
scalar[i+gidY*(WIDTH)] = ((distance(pos, center)<0.15f) ? 1.0f : 0.0f);
}
}
-__kernel void initVelocity(__global float* veloX,__global float* veloY,
- float4 minPos,
- float4 size)
+
+
+__kernel void velocity(__global float* veloX,__global float* veloY,
+ float t,
+ float4 minPos,
+ float4 size)
{
uint gidX = get_global_id(0);
uint gidY = get_global_id(1);
uint i;
float v;
+ float time_term = cos(t*M_PI_F/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);
+ 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;
veloX[i+gidY*(WIDTH)] = v; // = Vx(x,y)
veloY[i+gidY*(WIDTH)] = -v;// = Vy(y,x) = -Vx(x,y)
}
diff --git a/Examples/levelSet2D.py b/Examples/levelSet2D.py
index e5eb0d243..4fd1bcbd2 100644
--- a/Examples/levelSet2D.py
+++ b/Examples/levelSet2D.py
@@ -20,14 +20,14 @@ def scalaire(x, y):
def run():
dim = 2
- nb = 512
+ nb = 256
boxLength = [1., 1.]
boxMin = [0., 0.]
nbElem = [nb, nb]
timeStep = 0.05
#period = 10.
- finalTime = timeStep*2
+ finalTime = 3.
outputFilePrefix = './res/levelSet_2D_'
outputModulo = 1
@@ -44,12 +44,13 @@ def run():
## Operators
advec = pp.Transport(velo, scal)
+ 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, [advec])
+ pb = pp.Problem(topo3D, [velocity, advec])
## Setting solver to Problem
pb.setSolver(finalTime, timeStep, 'gpu',
diff --git a/Examples/levelSet3D.cl b/Examples/levelSet3D.cl
index 0f11c1a71..9a22b8ed8 100644
--- a/Examples/levelSet3D.cl
+++ b/Examples/levelSet3D.cl
@@ -15,30 +15,29 @@ __kernel void initScalar(__global float* scalar,
}
}
-__kernel void initVelocity(__global float* veloX,__global float* veloY,__global float* veloZ,
- float4 minPos,
- float4 size)
+__kernel void velocity(__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);
+
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;
- veloX[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = 2.0f * sin(pix)*sin(pix)*sin(2.0f*piy)*sin(2.0f*piz); // Vx(x,y,z)
- veloY[gidY+i*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0f*pix)*sin(piy)*sin(piy)*sin(2.0f*piz);// Vy(y,x,z)
- //veloY[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0f*piy)*sin(pix)*sin(pix)*sin(2.0f*piz);// Vy(y,x,z)
- veloZ[gidZ+i*(WIDTH)+gidY*(WIDTH*WIDTH)] = -sin(2.0f*pix)*sin(2.0f*piy)*sin(piz)*sin(piz);//Vz(z,x,y)
- //veloZ[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -sin(2.0f*piz)*sin(2.0f*pix)*sin(piy)*sin(piy);//Vz(z,x,y)
- /*
- v =-sin(piy)*sin(piy)*sin(2.0f*piz);
- veloX[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = 0.0f; // = Vx(x,y,z)
- veloY[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = v;// = Vy(y,x,z) = -Vx(x,y,z)
- veloZ[i+gidY*(WIDTH)+gidZ*(WIDTH*WIDTH)] = -v; // Vz(z,x,y) = -Vy(y,x,z)*/
+
+ 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)
}
}
diff --git a/Examples/levelSet3D.py b/Examples/levelSet3D.py
index 2db0da470..ecaf52bbd 100644
--- a/Examples/levelSet3D.py
+++ b/Examples/levelSet3D.py
@@ -13,7 +13,7 @@ def run():
timeStep = 0.05
#period = 10.
- finalTime = 20*timeStep#3.
+ finalTime = 3.
outputFilePrefix = './res/levelSet_3D_'
outputModulo = 1
@@ -30,12 +30,13 @@ def run():
## Operators
advec = pp.Transport(velo, scal)
+ 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, [advec])
+ pb = pp.Problem(topo3D, [velocity, advec])
## Setting solver to Problem
pb.setSolver(finalTime, timeStep, 'gpu',
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index 435cd62c9..dd25d5474 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -38,7 +38,9 @@ AnalyticalField = fields.analytical.AnalyticalField
## Operators
import operator.Transport
+import operator.Velocity
Transport = operator.Transport.Transport
+Velocity = operator.Velocity.Velocity
## Problem
import problem.problem
diff --git a/HySoP/hysop/operator/Velocity.py b/HySoP/hysop/operator/Velocity.py
new file mode 100644
index 000000000..f03384b5a
--- /dev/null
+++ b/HySoP/hysop/operator/Velocity.py
@@ -0,0 +1,52 @@
+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator.Advection
+
+Advection operator representation
+"""
+from .continuous import ContinuousOperator
+from .velocity_d import Velocity_P
+
+
+class Velocity(ContinuousOperator):
+ """
+ Velocity operator representation
+ """
+
+ def __init__(self, velocity):
+ """
+ Constructor.
+ Create an Advection operator from given variables velocity and scalar.
+
+ @param velocity ContinuousField : velocity variable.
+ @param scalar ContinuousField : scalar variable.
+ """
+ ContinuousOperator.__init__(self)
+ ## velocity variable
+ self.velocity = velocity
+ self.addVariable([velocity])
+
+ def discretize(self, idVelocityD=0, method=None):
+ """
+ Advection operator discretization method.
+ Create an AdvectionDOp.AdvectionDOp from given specifications.
+
+ @param spec : discretization specifications, not used.
+ """
+ if self.discreteOperator is None:
+ self.discreteOperator = Velocity_P(self, idVelocityD)
+
+ def __str__(self):
+ """ToString method"""
+ s = "Velocity operator (ContinuousOperator)"
+ if self.discreteOperator is not None:
+ s += " with the following discretization:\n"
+ s += str(self.discreteOperator)
+ else:
+ s += ". Not discretised"
+ return s + "\n"
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Advection"
+ print Advection.__doc__
diff --git a/HySoP/hysop/operator/advection_d.py b/HySoP/hysop/operator/advection_d.py
index 01a4b53c6..090c695ec 100644
--- a/HySoP/hysop/operator/advection_d.py
+++ b/HySoP/hysop/operator/advection_d.py
@@ -43,6 +43,7 @@ class Advection_P(DiscreteOperator):
self.compute_time = [0., 0., 0.]
self.compute_time_copy = [0., 0., 0.]
self.compute_time_transpose = [0., 0., 0.]
+ self.name = "advection"
def apply(self, t, dt, splittingDirection):
c_time, c_time_init = 0., 0.
diff --git a/HySoP/hysop/operator/remeshing.py b/HySoP/hysop/operator/remeshing.py
index 0bd0fa807..779aa93db 100644
--- a/HySoP/hysop/operator/remeshing.py
+++ b/HySoP/hysop/operator/remeshing.py
@@ -42,6 +42,7 @@ class Remeshing(DiscreteOperator):
self.output = [self.res_scalar]
self.needSplitting = True
self.compute_time = [0., 0., 0.]
+ self.name = "remeshing"
def apply(self, t, dt, splittingDirection):
"""
diff --git a/HySoP/hysop/operator/splitting.py b/HySoP/hysop/operator/splitting.py
index 286345303..c904e54d8 100644
--- a/HySoP/hysop/operator/splitting.py
+++ b/HySoP/hysop/operator/splitting.py
@@ -31,6 +31,7 @@ class Splitting(DiscreteOperator):
self.splitting.append((dim - 1, 1.))
[self.splitting.append((dim - 2 - i, 0.5)) for i in xrange(dim - 1)]
self.compute_time_details = [[0. for i in xrange(dim)] for op in self.operators]
+ self.name = "splitting"
def apply(self, t, dt):
"""
diff --git a/HySoP/hysop/operator/velocity_d.py b/HySoP/hysop/operator/velocity_d.py
new file mode 100644
index 000000000..2af16d6ae
--- /dev/null
+++ b/HySoP/hysop/operator/velocity_d.py
@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete velocity
+"""
+from ..constants import *
+from .discrete import DiscreteOperator
+import pyopencl as cl
+import numpy as np
+import time
+
+
+class Velocity_P(DiscreteOperator):
+ """
+ Velocity computation
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, advec, idVelocityD=0, method=None):
+ """
+ Constructor.
+ Create a Velocity operator on a given continuous domain.
+
+ """
+ DiscreteOperator.__init__(self)
+ ## Velocity.
+ self.velocity = advec.velocity.discreteField[idVelocityD]
+ ## input fields
+ self.input = [self.velocity]
+ ## output fields
+ self.output = [self.velocity]
+ self.numMethod = method
+ self.needSplitting = False
+ self.old_splitting_direction = None
+ self.compute_time = 0.
+ self.name = "velocity"
+
+ def apply(self, t, dt):
+ c_time = 0.
+ if self.numMethod is not None:
+ # Compute velocity field
+ dim = self.velocity.topology.dim
+ coord_min = np.ones(4, dtype=PARMES_REAL_GPU)
+ mesh_size = np.ones(4, dtype=PARMES_REAL_GPU)
+ coord_min[0:dim] = self.velocity.topology.mesh.start
+ mesh_size[0:dim] = self.velocity.topology.mesh.size
+ if len(self.velocity.gpu_data) == 3:
+ evt = self.numMethod.launch(self.velocity.gpu_data[0],
+ self.velocity.gpu_data[1],
+ self.velocity.gpu_data[2],
+ PARMES_REAL_GPU(t),
+ coord_min,
+ mesh_size)
+ else:
+ evt = self.numMethod.launch(self.velocity.gpu_data[0],
+ self.velocity.gpu_data[1],
+ PARMES_REAL_GPU(t),
+ coord_min,
+ mesh_size)
+ for df in self.output:
+ df.contains_data = False
+ self.numMethod.queue.finish()
+ c_time = (evt.profile.end - evt.profile.start) * 1e-9
+ self.compute_time += c_time
+ self.total_time += c_time
+ return c_time
+
+ def printComputeTime(self):
+ print "Velocity total time : ", self.total_time
+
+ def __str__(self):
+ s = "Velocity_P (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : AdvectionDOp"
+ print AdvectionDOp.__doc__
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index 7613d458b..4cdb74de8 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -80,7 +80,7 @@ class ParticleSolver(Solver):
else:
op.discretize()
## Create remeshing operator
- self.remesh = Remeshing(partPositions=self.advection.discreteOperator.res_position,
+ self.remeshing = Remeshing(partPositions=self.advection.discreteOperator.res_position,
partScalar=self.advection.discreteOperator.res_scalar,
resscal=self.advection.discreteOperator.scalar,
method=self.RemeshingMethod)
@@ -88,7 +88,7 @@ class ParticleSolver(Solver):
if op.discreteOperator.needSplitting:
if op is self.advection:
index = self.problem.operators.index(op)
- self.problem.operators[index] = Splitting([op, self.remesh], dim=self.problem.topology.dim)
+ self.problem.operators[index] = Splitting([op, self.remeshing], dim=self.problem.topology.dim)
else:
self.problem.operators[index] = Splitting([op], dim=self.problem.topology.dim)
self.isInitialized = True
diff --git a/HySoP/hysop/particular_solvers/gpu.py b/HySoP/hysop/particular_solvers/gpu.py
index fcf2dec8b..a8dd50c8a 100644
--- a/HySoP/hysop/particular_solvers/gpu.py
+++ b/HySoP/hysop/particular_solvers/gpu.py
@@ -223,7 +223,20 @@ class GPUParticleSolver(ParticleSolver):
if compute_time > 0.:
print "Total Computing : ", compute_time, "sec"
## Setting advection and remeshing kernels:
- self.advection.setMethod(KernelLauncher(self.prg.advection, self.queue, global_wg, local_wg))
+ for op in self.problem.operators:
+ if op.discreteOperator.name == "splitting":
+ for sop in op.operators:
+ print sop.discreteOperator.name
+ sop.setMethod(KernelLauncher(eval("self.prg."+sop.discreteOperator.name),
+ self.queue,
+ global_wg,
+ local_wg))
+ else:
+ print op.discreteOperator.name
+ op.setMethod(KernelLauncher(eval("self.prg."+op.discreteOperator.name),
+ self.queue,
+ global_wg,
+ local_wg))
self.advection.discreteOperator.init_copy = KernelLauncher(self.prg.advec_init_copy,
self.queue,
global_wg,
@@ -262,10 +275,6 @@ class GPUParticleSolver(ParticleSolver):
],
[(64, 1, 1),
(64, 1, 1)])
- self.remesh.setMethod(KernelLauncher(self.prg.remeshing,
- self.queue,
- global_wg,
- local_wg))
self.problem.io.set_get_data_method(self.get_data_from_device)
print "===\n"
--
GitLab