diff --git a/parmepy/new_ParMePy/Operator/AdvectionDOp.py b/parmepy/new_ParMePy/Operator/AdvectionDOp.py
index 1edac44dddac2e40cd3c5fc0947ab358374feb87..850f45469d6269d7e545274a304bce5a070528f0 100644
--- a/parmepy/new_ParMePy/Operator/AdvectionDOp.py
+++ b/parmepy/new_ParMePy/Operator/AdvectionDOp.py
@@ -61,9 +61,8 @@ class AdvectionDOp(DiscreteOperator):
c_time = time.time()
## ------------------- NumPy Optimisation
# Particle update
- self.resultScalar.values[:] = self.scalar.values[:]
- self.resultPosition.values[:] = self.velocity.domain.points[:]
- #print self.resultPosition.values
+ self.resultScalar.values[...] = self.scalar.values[...]
+ self.resultPosition.values[...] = self.velocity.domain.points[...]
# Advection
self.numMethod.integrate(self.velocity.domain.points, self.velocity.values, self.resultPosition.values, t, dt, splittingDirection)
#print self.resultPosition.values
diff --git a/parmepy/new_ParMePy/Operator/RemeshingDOp.py b/parmepy/new_ParMePy/Operator/RemeshingDOp.py
index 1a4e7d232f2a88d7bbda25f9d3f1031ccbfd8a6f..bb709ea38b16f9a378b92d9f91b818b8ff58235e 100644
--- a/parmepy/new_ParMePy/Operator/RemeshingDOp.py
+++ b/parmepy/new_ParMePy/Operator/RemeshingDOp.py
@@ -68,7 +68,7 @@ class RemeshingDOp(DiscreteOperator):
if self.gpu_kernel is None:
c_time = time.time()
## ------------------- NumPy Optimisation
- self.resultScalar.values[:] = 0.
+ self.resultScalar.values[...] = 0.
self.numMethod.remesh(self.ppos.values, self.pscal.values, splittingDirection)
## -------------------
self.compute_time += (time.time() - c_time)
diff --git a/parmepy/new_ParMePy/Problem/DiscreteTransportProblem.py b/parmepy/new_ParMePy/Problem/DiscreteTransportProblem.py
index 9b0d7f7681b2ff47ce37b1bf659284f2fdc3cef0..ddf8fbd93adbd01a9f71b101e2c1286a78dbaee7 100644
--- a/parmepy/new_ParMePy/Problem/DiscreteTransportProblem.py
+++ b/parmepy/new_ParMePy/Problem/DiscreteTransportProblem.py
@@ -66,10 +66,10 @@ class DiscreteTransportProblem(DiscreteProblem):
for op in self.operators:
c_time += op.compute_time
c_flop += op.total_flop
- print "Time", op, op.compute_time, "s \t", (op.total_flop/op.compute_time) * 1e-9, "GFlop/s"
+ print "Time", op, op.compute_time, "s \t", (op.total_flop / op.compute_time) * 1e-9, "GFlop/s"
q_time += op.queued_time
s_time += op.submit_time
- print "Computing time : ", c_time, "s \t", (c_flop/c_time) * 1e-9, "GFlop/s"
+ print "Computing time : ", c_time, "s \t", (c_flop / c_time) * 1e-9, "GFlop/s"
if not self.operators[0].gpu_kernel is None:
print "OpenCL queued time : ", q_time, "s"
print "OpenCL submit time : ", s_time, "s"
diff --git a/parmepy/new_ParMePy/Utils/Linear.py b/parmepy/new_ParMePy/Utils/Linear.py
index be85e2eb62bc8064423d2c0377236ad9382f42c1..38fe792b5c427f412eb05b5c657bcb078f1442b2 100644
--- a/parmepy/new_ParMePy/Utils/Linear.py
+++ b/parmepy/new_ParMePy/Utils/Linear.py
@@ -25,7 +25,6 @@ class Linear(InterpolationMethod):
self.grid = grid
self.gvalues = gvalues
self.dim_range = xrange(self.grid.dimension)
- self.ind = np.empty(grid.shape, dtype=dtype_integer)
self.nb_flop = 9
def interpolate(self, t, ppos, dir):
@@ -38,9 +37,12 @@ class Linear(InterpolationMethod):
@param dir : interpolation direction
@return particle velocity, velocity at position ppos.
"""
+ ind = np.empty(self.grid.shape, dtype=dtype_integer)
a0 = ((ppos - self.grid.min) / self.grid.elementSize)
- self.ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
- index = [self.ind[..., i] for i in self.dim_range]
+ ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
+ ind[..., dir] = np.floor(a0[..., dir])
+ #ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
+ index = [ind[..., i] for i in self.dim_range]
index[dir] = (index[dir] + self.grid.elementNumber[dir]) % self.grid.elementNumber[dir]
y = (ppos[..., dir] - self.grid[index][..., dir]) / self.grid.elementSize[dir]
res = np.copy(self.gvalues.values)
diff --git a/parmepy/new_ParMePy/Utils/M6Prime.py b/parmepy/new_ParMePy/Utils/M6Prime.py
index 98f6f289b7b6e62eb10c4ec01ae733ba36b615e0..e965760bc93e84c42be8e82335a523bc01549df4 100644
--- a/parmepy/new_ParMePy/Utils/M6Prime.py
+++ b/parmepy/new_ParMePy/Utils/M6Prime.py
@@ -17,7 +17,6 @@ class M6Prime(RemeshingMethod):
"""
RemeshingMethod.__init__(self, grid, gvalues)
self.dim_range = xrange(self.grid.dimension)
- self.ind = np.empty(grid.shape, dtype=dtype_integer)
self.nb_flop = 106
def remesh(self, ppos, pscal, dir):
@@ -28,28 +27,41 @@ class M6Prime(RemeshingMethod):
@param pscal : particle scalar.
@param dir : direction to remesh.
"""
+ ind = np.empty(self.grid.shape, dtype=dtype_integer)
a0 = ((ppos - self.grid.min) / self.grid.elementSize)
- self.ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
- index = [self.ind[..., i] for i in self.dim_range]
- x0 = (ppos[..., dir] - self.ind[..., dir] * self.grid.elementSize[dir] - self.grid.min[dir]) / self.grid.elementSize[dir]
+ ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
+ ind[..., dir] = np.floor(a0[..., dir])
+ x0 = (ppos[..., dir] - ind[..., dir] * self.grid.elementSize[dir] - self.grid.min[dir]) / self.grid.elementSize[dir]
am2 = (-(x0) * (5. * (x0 + 2.) - 8.) * (x0 - 1.) ** 3) / 24.
am = x0 * (x0 - 1.) * (25. * (x0 + 1.) ** 3 - 114. * (x0 + 1.) ** 2 + 153. * (x0 + 1.) - 48.) / 24.
a0 = -(x0 - 1.) * (25. * x0 ** 4 - 38. * x0 ** 3 - 3. * x0 ** 2 + 12. * x0 + 12.) / 12.
ap = x0 * (25. * (1. - x0) ** 4 - 38. * (1. - x0) ** 3 - 3. * (1. - x0) ** 2 + 12. * (1. - x0) + 12.) / 12.
ap2 = (1. - x0) * (-x0) * (25. * (2. - x0) ** 3 - 114. * (2. - x0) ** 2 + 153. * (2. - x0) - 48.) / 24.
ap3 = -(1. - x0) * ((5. * (3. - x0) - 8.) * (-x0) ** 3) / 24.
- index[dir] = (index[dir] - 2 + self.grid.elementNumber[dir]) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * am2
- index[dir] = (index[dir] + 1 ) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * am
- index[dir] = (index[dir] + 1 ) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * a0
- index[dir] = (index[dir] + 1 ) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * ap
- index[dir] = (index[dir] + 1 ) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * ap2
- index[dir] = (index[dir] + 1 ) % self.grid.elementNumber[dir]
- self.gvalues.values[tuple(index)] += pscal * ap3
+
+ ind[..., dir] = (ind[..., dir] - 2 + self.grid.elementNumber[dir]) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), am2.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
+
+ ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), am.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
+
+ ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), a0.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
+
+ ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
+
+ ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap2.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
+
+ ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
+ for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap3.reshape((-1))):
+ self.gvalues.values[tuple(i_index)] += ps * w
def __str__(self):
"""ToString method"""