diff --git a/HySoP/hysop/fields/scalar.py b/HySoP/hysop/fields/scalar.py
index 02f69b4f3d6aa4e79e8c69858e0c94f94169a099..0140004776d7451ea3714e6519ca0d06fbc0a051 100644
--- a/HySoP/hysop/fields/scalar.py
+++ b/HySoP/hysop/fields/scalar.py
@@ -82,7 +82,7 @@ class ScalarField(DiscreteField):
else:
self.data[...] = np.asarray(
- formula(*(self.topology.mesh.coords + args)),
+ v_formula(self.topology.mesh.coords, *args),
dtype=PARMES_REAL,
order=ORDER)
else:
diff --git a/HySoP/hysop/gpu/gpu_particle_advection.py b/HySoP/hysop/gpu/gpu_particle_advection.py
index 53d06030f3a71870c2265a1598704721c70a3001..8d23b5798dfe5c0f0471d677816c768bfc5e07b6 100644
--- a/HySoP/hysop/gpu/gpu_particle_advection.py
+++ b/HySoP/hysop/gpu/gpu_particle_advection.py
@@ -10,8 +10,8 @@ from parmepy.operator.discrete.particle_advection import ParticleAdvection
from parmepy.gpu import PARMES_REAL_GPU, cl
from parmepy.gpu.tools import get_opencl_environment
from parmepy.gpu.gpu_kernel import KernelLauncher
-from parmepy.gpu.gpu_kernel import KernelListLauncher
-from parmepy.numerics.splitting import Splitting
+#from parmepy.gpu.gpu_kernel import KernelListLauncher
+#from parmepy.numerics.splitting import Splitting
class GPUParticleAdvection(ParticleAdvection):
diff --git a/HySoP/hysop/mpi/mesh.py b/HySoP/hysop/mpi/mesh.py
index 99b8fdc1c9dcdb69db20557a3915577c9e29913a..9137455380ba3f9e26cfe2a150cc71c9850c252a 100644
--- a/HySoP/hysop/mpi/mesh.py
+++ b/HySoP/hysop/mpi/mesh.py
@@ -56,7 +56,8 @@ class SubMesh(object):
self.end = self.origin + self.space_step * (self.resolution - 1)
if self.dim == 1:
- self.coords = np.linspace(self.origin, self.end, self.resolution)
+ cx = np.linspace(self.origin, self.end, self.resolution)
+ self.coords = tuple([cx,])
elif self.dim == 2:
cx = np.linspace(self.origin[0], self.end[0],
self.resolution[0])[:, np.newaxis]
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta2.py b/HySoP/hysop/numerics/integrators/runge_kutta2.py
index 24e6bec3592d9cff118a4f554ada5654b4794bae..80471ce9f434d2efbedea22836f35469e45139be 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta2.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta2.py
@@ -74,7 +74,8 @@ class RK2_scalar(ODESolver):
ODESolver.__init__(self, f)
self.name = 'RK2'
- def __call__(self, f, t, dt, y, res=None, work=None, **f_kwargs):
+ def __call__(self, f, t, dt, y, res=None, work=None, yprime=None,
+ **f_kwargs):
"""
Integration step for RK2 Method.
@@ -86,13 +87,17 @@ class RK2_scalar(ODESolver):
@param y : function Y(t,space).
@param res : result (may be allocated when not None).
@param work : allocated working array.
+ @param yprime : y'(t) already computed (avoid 1st call to f)
@param f_kwargs : Parameters for right-hand side computing.
"""
if res is None:
res = np.empty_like(y)
if work is None:
work = np.empty_like(y)
- work[...] = f(t, y, res=work, **f_kwargs)
+ if yprime is None:
+ work[...] = f(t, y, res=work, **f_kwargs)
+ else:
+ work[...] = yprime
work[...] *= 0.5 * dt
res[...] = y
res[...] += dt * f(t + 0.5 * dt, y + work, **f_kwargs)
diff --git a/HySoP/hysop/numerics/interpolation.py b/HySoP/hysop/numerics/interpolation.py
index 3561ac2ae0ef18472345c384648e0a1673391f8f..246e72f54336b586e0a3c9ae3c09fc169ab6d386 100644
--- a/HySoP/hysop/numerics/interpolation.py
+++ b/HySoP/hysop/numerics/interpolation.py
@@ -39,9 +39,10 @@ class Linear(NumMethod):
i_row, i_col = np.indices(x.shape)
if self.d == 0:
floor = res
- floor[...] = np.floor(
- (x - self.origin[self.d]) / self.dx[self.d])
- i_y[...] = (x - self.origin[self.d]) / self.dx[self.d] - floor
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
i_row[...] = floor.astype(np.int32) % x.shape[self.d]
i_col[...] = np.indices((x.shape[1],))[0][np.newaxis, :]
res[...] = self.field.data[self.d][i_row, i_col] * (1. - i_y)
@@ -49,9 +50,10 @@ class Linear(NumMethod):
res[...] += self.field.data[self.d][i_row, i_col] * i_y
else:
floor = res
- floor[...] = np.floor(
- (x - self.origin[self.d]) / self.dx[self.d])
- i_y[...] = (x - self.origin[self.d]) / self.dx[self.d] - floor
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
i_col[...] = floor.astype(np.int32) % x.shape[self.d]
i_row[...] = np.indices((x.shape[0],))[0][:, np.newaxis]
res[...] = self.field.data[self.d][i_row, i_col] * (1. - i_y)
@@ -63,9 +65,10 @@ class Linear(NumMethod):
i_row, i_col, i_dep = np.indices(x.shape)
if self.d == 0:
floor = res
- floor[...] = np.floor(
- (x - self.origin[self.d]) / self.dx[self.d])
- i_y[...] = (x - self.origin[self.d]) / self.dx[self.d] - floor
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
i_row[...] = floor.astype(np.int32) % x.shape[self.d]
i_col[...] = \
np.indices((x.shape[1],))[0][np.newaxis, :, np.newaxis]
@@ -78,9 +81,10 @@ class Linear(NumMethod):
self.field.data[self.d][i_row, i_col, i_dep] * i_y
elif self.d == 1:
floor = res
- floor[...] = np.floor(
- (x - self.origin[self.d]) / self.dx[self.d])
- i_y[...] = (x - self.origin[self.d]) / self.dx[self.d] - floor
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
i_col[...] = floor.astype(np.int32) % x.shape[self.d]
i_row[...] = \
np.indices((x.shape[0],))[0][:, np.newaxis, np.newaxis]
@@ -93,9 +97,10 @@ class Linear(NumMethod):
self.field.data[self.d][i_row, i_col, i_dep] * i_y
else:
floor = res
- floor[...] = np.floor(
- (x - self.origin[self.d]) / self.dx[self.d])
- i_y[...] = (x - self.origin[self.d]) / self.dx[self.d] - floor
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
i_dep[...] = floor.astype(np.int32) % x.shape[self.d]
i_row[...] = \
np.indices((x.shape[0],))[0][:, np.newaxis, np.newaxis]
@@ -107,4 +112,17 @@ class Linear(NumMethod):
res[...] += \
self.field.data[self.d][i_row, i_col, i_dep] * i_y
+ else:
+ if i_row is None:
+ i_row = np.indices(x.shape)
+ floor = res
+ floor[...] = (x - self.origin[self.d]) / self.dx[self.d]
+ i_y[...] = floor
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
+ i_row[...] = floor.astype(np.int32) % x.shape[self.d]
+ res[...] = self.field.data[i_row] * (1. - i_y)
+ i_row[...] = (i_row + 1) % x.shape[self.d]
+ res[...] += self.field.data[i_row] * i_y
+
return res
diff --git a/HySoP/hysop/numerics/remeshing/l6star.py b/HySoP/hysop/numerics/remeshing/l6star.py
deleted file mode 100644
index 6f179ea4edb2f3b9adbedba721ca82ffafd52c9c..0000000000000000000000000000000000000000
--- a/HySoP/hysop/numerics/remeshing/l6star.py
+++ /dev/null
@@ -1,187 +0,0 @@
-"""
-@file l6star.py
-"""
-from parmepy.constants import np, PARMES_REAL, ORDER
-from parmepy.numerics.remeshing.remeshing import Remeshing
-
-
-class L6Star(Remeshing):
- """Remshing with Lambda6star function"""
-
- def __init__(self, dim, dx, origin):
- Remeshing.__init__(self, dim)
- self.name = 'Lambda6star'
- self.dx = dx
- self.origin = origin
- self._alpha = lambda y, s: s * (((-12. + (4. + (15. + (140. + (-370.
- + (312. - 89. * y) * y) * y) * y) * y) * y) * y) / 720.)
- self._beta = lambda y, s: s * (((108. + (-54. + (-120. + (-955. +
- (2581. + (-2183. + 623. * y) * y) * y) * y) * y) * y ) * y) / 720.)
- self._gamma = lambda y, s: s * (((-180. + (180. + (65. + (950. +
- (-2574. + (2182. - 623. * y) * y) * y) * y) * y) * y) * y) / 240.)
- self._delta = lambda y, s: s * ( 1. + ((-196. + (-959. + (2569. +
- (-2181. + 623. * y) * y) * y) * y * y ) * y * y) / 144.)
- self._eta = lambda y, s: s * (((108. + (108. + (-39. + (976. + (-2566.
- + (2180. - 623. * y) * y) * y) * y) * y) * y) * y) / 144.)
- self._zeta = lambda y, s: s * (((-36. + (-18. + (40. + (-995. + (2565.
- + (-2179. + 623. * y) * y) * y) * y) * y) * y) * y) / 240.)
- self._theta = lambda y, s: s * (((12. + (4. + (-15. + (1010. + (-2566.
- + (2178. - 623. * y) * y) * y) * y) * y) * y) * y) / 720.)
- self._iota = lambda y, s: s * (((-145. + (367. + (-311. + 89. * y)
- * y) * y) * y * y * y * y) / 720.)
-
- def __call__(self, ppos, pscal, d):
-
- """
- Remesh particles at position ppos with scalar
- pscal along direction d.
-
- @param ppos : particle position
- @param pscal : particle scalar
- @param d : splitting direction
- """
- if len(ppos.shape) == 2:
- size = ppos.shape
- index = np.indices(size)
- floor = np.floor(ppos / self.dx[d])
- y = np.array(ppos / self.dx[d] - floor,
- dtype=PARMES_REAL, order=ORDER)
- result = np.zeros(size, dtype=PARMES_REAL, order=ORDER)
- ind = floor.astype(np.int32) % size[d]
-
- # alpha
- index[d] = (ind - 3) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._alpha(y, pscal)[sl]
-
- # beta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._beta(y, pscal)[sl]
-
- # gamma
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._gamma(y, pscal)[sl]
-
- # delta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._delta(y, pscal)[sl]
-
- # eta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._eta(y, pscal)[sl]
-
- # zeta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._zeta(y, pscal)[sl]
-
- # theta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._theta(y, pscal)[sl]
-
- # iota
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl]] = \
- result[index[0][sl], index[1][sl]] +\
- self._iota(y, pscal)[sl]
- else:
- size = ppos.shape
- index = np.indices(size)
- floor = np.floor(ppos / self.dx[d])
- y = np.array(ppos / self.dx[d] - floor,
- dtype=PARMES_REAL, order=ORDER)
- result = np.zeros(size, dtype=PARMES_REAL, order=ORDER)
- ind = floor.astype(np.int32) % size[d]
-
- # alpha
- index[d] = (ind - 3) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] + \
- self._alpha(y, pscal)[sl]
-
- # beta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._beta(y, pscal)[sl]
-
- # gamma
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._gamma(y, pscal)[sl]
-
- # delta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] + \
- self._delta(y, pscal)[sl]
-
- # eta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._eta(y, pscal)[sl]
-
- # zeta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._zeta(y, pscal)[sl]
-
- # theta
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._theta(y, pscal)[sl]
-
- # iota
- index[d] = (index[d] + 1) % size[d]
- for i in xrange(size[d]):
- sl = self.slice_i_along_d(i, d)
- result[index[0][sl], index[1][sl], index[2][sl]] = \
- result[index[0][sl], index[1][sl], index[2][sl]] +\
- self._iota(y, pscal)[sl]
- return result
diff --git a/HySoP/hysop/numerics/remeshing/m4prime.py b/HySoP/hysop/numerics/remeshing/m4prime.py
deleted file mode 100644
index 949edca08f1d5ffec6092090c975d9492c498f54..0000000000000000000000000000000000000000
--- a/HySoP/hysop/numerics/remeshing/m4prime.py
+++ /dev/null
@@ -1,171 +0,0 @@
-"""
-@file m4prime.py
-"""
-from parmepy.constants import np, PARMES_REAL, ORDER
-from parmepy.numerics.remeshing.remeshing import Remeshing
-
-
-class M4Prime(Remeshing):
- """Remshing with M4prime function"""
-
- def __init__(self, dim, dx, origin):
- Remeshing.__init__(self, dim)
- self.name = 'M4Prime'
- self.dx = dx
- self.origin = origin
- self._alpha = lambda y, s: s * ((y * (y * (-y + 2.) - 1.)) / 2.)
- self._beta = lambda y, s: s * ((y * y * (3. * y - 5.) + 2.) / 2.)
- self._gamma = lambda y, s: s * ((y * (y * (-3. * y + 4.) + 1.)) / 2.)
- self._delta = lambda y, s: s * ((y * y * (y - 1.)) / 2.)
-
- def __call__(self, ppos, origin, pscal, d, res=None,
- i_row=None, i_col=None, i_dep=None, i_y=None):
-
- """
- Remesh particles at position ppos with scalar
- pscal along direction d.
-
- @param ppos : particle position
- @param pscal : particle scalar
- @param d : splitting direction
- @param i_row : allocated working array (rows indices array).
- @param i_col : allocated working array (columns indices array).
- @param i_dep : allocated working array (depth indices array).
- @param i_y : allocated working array (distance to left grid point).
- """
- shape = ppos.shape
- if len(ppos.shape) == 2:
- index = [i_row, i_col]
- if d == 0:
- index[1][...] = np.indices((shape[1],))[0][np.newaxis, :]
- else:
- index[0][...] = np.indices((shape[0],))[0][:, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 1) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._delta(i_y, pscal)[sl]
-
- elif len(ppos.shape) == 3:
- index = [i_row, i_col, i_dep]
- if d == 0:
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- elif d == 1:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- else:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 1) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._delta(i_y, pscal)[sl]
-
- else:
- index = [i_row]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 1) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._delta(i_y, pscal)[sl]
diff --git a/HySoP/hysop/numerics/remeshing/m6prime.py b/HySoP/hysop/numerics/remeshing/m6prime.py
deleted file mode 100644
index e320b373313fb0267cd51ea1a81387490fc7b2db..0000000000000000000000000000000000000000
--- a/HySoP/hysop/numerics/remeshing/m6prime.py
+++ /dev/null
@@ -1,225 +0,0 @@
-"""
-@file m6prime.py
-"""
-from parmepy.constants import np
-from parmepy.numerics.remeshing.remeshing import Remeshing
-
-
-class M6Prime(Remeshing):
- """Remshing with M6prime function"""
-
- def __init__(self, dim, dx, origin):
- Remeshing.__init__(self, dim)
- self.name = 'M6Prime'
- self.dx = dx
- self.origin = origin
- self._alpha = lambda y, s: s * y * (
- y * (y * (y * (13. - 5. * y) - 9.) - 1.) + 2.) / 24.
- self._beta = lambda y, s: s * y * (
- y * (y * (y * (25. * y - 64.) + 39.) + 16.) - 16.) / 24.
- self._gamma = lambda y, s: s * (
- y * y * (y * (y * (126. - 50. * y) - 70.) - 30.) / 24. + 1.)
- self._delta = lambda y, s: s * y * (
- y * (y * (y * (50. * y - 124.) + 66.) + 16.) + 16.) / 24.
- self._eta = lambda y, s: s * y * (
- y * (y * (y * (61. - 25. * y) - 33.) - 1.) - 2.) / 24.
- self._zeta = lambda y, s: s * y * y * y * (
- y * (5. * y - 12.) + 7.) / 24.
-
- def __call__(self, ppos, origin, pscal, d, res=None,
- i_row=None, i_col=None, i_dep=None, i_y=None):
-
- """
- Remesh particles at position ppos with scalar
- pscal along direction d.
-
- @param ppos : particle position
- @param pscal : particle scalar
- @param d : splitting direction
- @param i_row : allocated working array (rows indices array).
- @param i_col : allocated working array (columns indices array).
- @param i_dep : allocated working array (depth indices array).
- @param i_y : allocated working array (distance to left grid point).
- """
- shape = ppos.shape
- if len(shape) == 2:
- index = [i_row, i_col]
- if d == 0:
- index[1][...] = np.indices((shape[1],))[0][np.newaxis, :]
- else:
- index[0][...] = np.indices((shape[0],))[0][:, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 2) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._zeta(i_y, pscal)[sl]
- elif len(ppos.shape) == 3:
- index = [i_row, i_col, i_dep]
- if d == 0:
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- elif d == 1:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- else:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 2) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._zeta(i_y, pscal)[sl]
- else:
- index = [i_row]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 2) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._zeta(i_y, pscal)[sl]
diff --git a/HySoP/hysop/numerics/remeshing/m8prime.py b/HySoP/hysop/numerics/remeshing/m8prime.py
deleted file mode 100644
index e5527e31a025064a2ca3d7644250b8d9237b61e5..0000000000000000000000000000000000000000
--- a/HySoP/hysop/numerics/remeshing/m8prime.py
+++ /dev/null
@@ -1,286 +0,0 @@
-"""
-@file m8prime.py
-"""
-from parmepy.constants import np
-from parmepy.numerics.remeshing.remeshing import Remeshing
-
-
-class M8Prime(Remeshing):
- """Remshing with M8prime function"""
-
- def __init__(self, dim, dx, origin):
- Remeshing.__init__(self, dim)
- self.name = 'M8Prime'
- self.dx = dx
- self.origin = origin
- self._alpha = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (-10. * y + 21.) + 28.) - 105.)
- + 70.) + 35.) - 56.) + 17.) / 3360.)
- self._beta = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (70. * y - 175.) - 140.) + 770.)
- - 560.) - 350.) + 504.) - 102.) / 3360.)
- self._gamma = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (-210. * y + 609.) + 224.) - 2135.)
- + 910.) + 2765.) - 2520.) + 255.) / 3360.)
- self._delta = lambda y, s: s * (
- (y * y * (y * y * (y * y * (70. * y - 231.) + 588.) - 980.) + 604.)
- / 672.)
- self._eta = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (-70. * y + 259.) - 84.) - 427.)
- - 182.) + 553.) + 504.) + 51.) / 672.)
- self._zeta = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (210. * y - 861.) + 532.) + 770.)
- + 560.) - 350.) - 504.) - 102.) / 3360.)
- self._theta = lambda y, s: s * (
- (y * (y * (y * (y * (y * (y * (-70. * y + 315.) - 280.) - 105.)
- - 70.) + 35.) + 56.) + 17.) / 3360.)
- self._iota = lambda y, s: s * (
- (y * y * y * y * y * (y * (10. * y - 49.) + 56.)) / 3360.)
-
- def __call__(self, ppos, origin, pscal, d, res=None,
- i_row=None, i_col=None, i_dep=None, i_y=None):
-
- """
- Remesh particles at position ppos with scalar
- pscal along direction d.
-
- @param ppos : particle position
- @param pscal : particle scalar
- @param d : splitting direction
- @param i_row : allocated working array (rows indices array).
- @param i_col : allocated working array (columns indices array).
- @param i_dep : allocated working array (depth indices array).
- @param i_y : allocated working array (distance to left grid point).
- """
- shape = ppos.shape
- if len(shape) == 2:
- index = [i_row, i_col]
- if d == 0:
- index[1][...] = np.indices((shape[1],))[0][np.newaxis, :]
- else:
- index[0][...] = np.indices((shape[0],))[0][:, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 3) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._zeta(i_y, pscal)[sl]
-
- # theta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._theta(i_y, pscal)[sl]
-
- # iota
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl]] = \
- res[index[0][sl], index[1][sl]] +\
- self._iota(i_y, pscal)[sl]
-
- elif len(ppos.shape) == 3:
- index = [i_row, i_col, i_dep]
- if d == 0:
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- elif d == 1:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[2][...] = \
- np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
- else:
- index[0][...] = \
- np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
- index[1][...] = \
- np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 3) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] + \
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._zeta(i_y, pscal)[sl]
-
- # theta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._theta(i_y, pscal)[sl]
-
- # iota
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl], index[1][sl], index[2][sl]] = \
- res[index[0][sl], index[1][sl], index[2][sl]] +\
- self._iota(i_y, pscal)[sl]
-
- else:
- index = [i_row]
- floor = res
- floor[...] = np.floor((ppos - origin) / self.dx[d])
- i_y[...] = ((ppos - origin) / self.dx[d]) - floor
-
- # alpha
- index[d][...] = (floor.astype(np.int32) - 3) % shape[d]
- res[...] = 0.
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._alpha(i_y, pscal)[sl]
-
- # beta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._beta(i_y, pscal)[sl]
-
- # gamma
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._gamma(i_y, pscal)[sl]
-
- # delta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._delta(i_y, pscal)[sl]
-
- # eta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._eta(i_y, pscal)[sl]
-
- # zeta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._zeta(i_y, pscal)[sl]
-
- # theta
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._theta(i_y, pscal)[sl]
-
- # iota
- index[d][...] = (index[d] + 1) % shape[d]
- for i in xrange(shape[d]):
- sl = self.slice_i_along_d(i, d)
- res[index[0][sl]] = \
- res[index[0][sl]] +\
- self._iota(i_y, pscal)[sl]
diff --git a/HySoP/hysop/numerics/remeshing/remeshing.py b/HySoP/hysop/numerics/remeshing/remeshing.py
index 8ef86e1c9adc5bbffdad65b90efc6ef65fc12ff8..ec05f98c21392ab048ae15b35163d275b0205458 100644
--- a/HySoP/hysop/numerics/remeshing/remeshing.py
+++ b/HySoP/hysop/numerics/remeshing/remeshing.py
@@ -1,30 +1,142 @@
"""
@file remeshing.py
"""
-from abc import ABCMeta, abstractmethod
-from parmepy.constants import np, PARMES_REAL, ORDER
+from parmepy.constants import np
from parmepy.numerics.method import NumMethod
class Remeshing(NumMethod):
- """Remshing with M6prime function"""
+ """Remshing"""
- __metaclass__ = ABCMeta
-
- @abstractmethod
- def __init__(self, dim):
+ def __init__(self, dim, dx, origin, formula):
+ """
+ Create a remeshing numeric method based on given formula.
+ @param dim : problem dimension
+ @param dx : mesh space step
+ @param origin : mesh lower point
+ @param formula : Remeshing formula to use.
+ Availables formulas :
+ - M4prime
+ - M6prime
+ - M8prime
+ - L6star
+ """
NumMethod.__init__(self)
- self.name = 'M6Prime'
+ self.name = formula
+ self.dx = dx
+ self.origin = origin
self._slice_all = [slice(None, None, None)
for d in xrange(dim)]
+ if self.name == 'm4prime':
+ self._shift = 1
+ self._weights = [
+ lambda y, s: s * ((y * (y * (-y + 2.) - 1.)) / 2.),
+ lambda y, s: s * ((y * y * (3. * y - 5.) + 2.) / 2.),
+ lambda y, s: s * ((y * (y * (-3. * y + 4.) + 1.)) / 2.),
+ lambda y, s: s * ((y * y * (y - 1.)) / 2.)
+ ]
+ elif self.name == 'l4star':
+ self._shift = 2
+ self._weights = [
+ lambda y, s: s * y * (y * (y * (y * (13. - 5. * y) - 9.) - 1.) + 2.) / 24.,
+ lambda y, s: s * y * (y * (y * (y * (25. * y - 64.) + 39.) + 16.) - 16.) / 24.,
+ lambda y, s: s * (y * y * (y * (y * (126. - 50. * y) - 70.) - 30.) / 24. + 1.),
+ lambda y, s: s * y * (y * (y * (y * (50. * y - 124.) + 66.) + 16.) + 16.) / 24.,
+ lambda y, s: s * y * (y * (y * (y * (61. - 25. * y) - 33.) - 1.) - 2.) / 24.,
+ lambda y, s: s * y * y * y * (y * (5. * y - 12.) + 7.) / 24.
+ ]
+ elif self.name == 'l4star_c3':
+ self._shift = 2
+ self._weights = [
+ lambda y, s: s * (2. + (-1. + (-2. + (-22. + (58. + (-49. + 14. * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (-16. + (16. + (4. + (111. + (-290. + (245. - 70. * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (1. + (-30. + (-224. + (580. + (-490. + 140. * y) * y) * y) * y * y) * y * y / 24.),
+ lambda y, s: s * (16. + (16. + (-4. + (226. + (-580. + (490. - 140. * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (-2. + (-1. + (2. + (-114. + (290. + (-245. + 70. * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (23. + (-58. + (49. - 14. * y) * y) * y) * y * y * y * y / 24.
+
+ ]
+ elif self.name == 'l4star_c4':
+ self._shift = 2
+ self._weights = [
+ lambda y, s: s * (2. + (-1. + (-2. + (1. + (-80. + (273. + (-354. + (207. - 46. * y) * y) * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (-16. + (16. + (4. + (-4. + (400. + (-1365. + (1770. + (-1035. + 230. * y) * y) * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (1. + (-30. + (6. + (-800. + (2730. + (-3540. + (2070. - 460. * y) * y) * y) * y) * y) * y * y) * y * y / 24.),
+ lambda y, s: s * (16. + (16. + (-4. + (-4. + (800. + (-2730. + (3540. + (-2070. + 460. * y) * y) * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (-2. + (-1. + (2. + (1. + (-400. + (1365. + (-1770. + (1035. - 230. * y) * y) * y) * y) * y) * y) * y) * y) * y / 24.,
+ lambda y, s: s * (80. + (-273. + (354. + (-207. + 46. * y) * y) * y) * y) * y * y * y * y * y / 24.
+ ]
+ elif self.name == 'm8prime':
+ self._shift = 3
+ self._weights = [
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (-10. * y + 21.) + 28.) - 105.)+ 70.) + 35.) - 56.) + 17.) / 3360.),
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (70. * y - 175.) - 140.) + 770.)- 560.) - 350.) + 504.) - 102.) / 3360.),
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (-210. * y + 609.) + 224.) - 2135.)+ 910.) + 2765.) - 2520.) + 255.) / 3360.),
+ lambda y, s: s * ((y * y * (y * y * (y * y * (70. * y - 231.) + 588.) - 980.) + 604.)/ 672.),
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (-70. * y + 259.) - 84.) - 427.)- 182.) + 553.) + 504.) + 51.) / 672.),
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (210. * y - 861.) + 532.) + 770.)+ 560.) - 350.) - 504.) - 102.) / 3360.),
+ lambda y, s: s * ((y * (y * (y * (y * (y * (y * (-70. * y + 315.) - 280.) - 105.)- 70.) + 35.) + 56.) + 17.) / 3360.),
+ lambda y, s: s * ((y * y * y * y * y * (y * (10. * y - 49.) + 56.)) / 3360.)
+ ]
+ elif self.name == 'l6star':
+ self._shift = 3
+ self._weights = [
+ lambda y, s: s * (((-12. + (4. + (15. + (140. + (-370.+ (312. - 89. * y) * y) * y) * y) * y) * y) * y) / 720.),
+ lambda y, s: s * (((108. + (-54. + (-120. + (-955. +(2581. + (-2183. + 623. * y) * y) * y) * y) * y) * y ) * y) / 720.),
+ lambda y, s: s * (((-180. + (180. + (65. + (950. +(-2574. + (2182. - 623. * y) * y) * y) * y) * y) * y) * y) / 240.),
+ lambda y, s: s * ( 1. + ((-196. + (-959. + (2569. +(-2181. + 623. * y) * y) * y) * y * y ) * y * y) / 144.),
+ lambda y, s: s * (((108. + (108. + (-39. + (976. + (-2566.+ (2180. - 623. * y) * y) * y) * y) * y) * y) * y) / 144.),
+ lambda y, s: s * (((-36. + (-18. + (40. + (-995. + (2565.+ (-2179. + 623. * y) * y) * y) * y) * y) * y) * y) / 240.),
+ lambda y, s: s * (((12. + (4. + (-15. + (1010. + (-2566.+ (2178. - 623. * y) * y) * y) * y) * y) * y) * y) / 720.),
+ lambda y, s: s * (((-145. + (367. + (-311. + 89. * y)* y) * y) * y * y * y * y) / 720.)
+ ]
+ elif self.name == 'l6star_c4':
+ self._shift = 3
+ self._weights = [
+ lambda y, s: s * (-12. + (4. + (15. + (-5. + (500. + (-1718. + (2231. + (-1305. + 290. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (108. + (-54. + (-120. + (60. + (-3509. + (12027. + (-15617. + (9135. - 2030. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-540. + (540. + (195. + (-195. + (10548. + (-36084. + (46851. + (-27405. + 6090. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (1. + (-980. + (280. + (-17605. + (60145. + (-78085. + (45675. - 10150. * y) * y) * y) * y) * y) * y * y) * y * y / 720.),
+ lambda y, s: s * (540. + (540. + (-195. + (-195. + (17620. + (-60150. + (78085. + (-45675. + 10150. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-108. + (-54. + (120. + (60. + (-10575. + (36093. + (-46851. + (27405. - 6090. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (12. + (4. + (-15. + (-5. + (3524. + (-12032. + (15617. + (-9135. + 2030. * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-503. + (1719. + (-2231. + (1305. - 290. * y) * y) * y) * y) * y * y * y * y * y / 720.
+ ]
+ elif self.name == 'l6star_c5':
+ self._shift = 3
+ self._weights = [
+ lambda y, s: s * (-12. + (4. + (15. + (-5. + (-3. + (1803. + (-7829. + (13785. + (-12285. + (5533. - 1006. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (108. + (-54. + (-120. + (60. + (12. + (-12620. + (54803. + (-96495. + (85995. + (-38731. + 7042. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-540. + (540. + (195. + (-195. + (-15. + (37857. + (-164409. + (289485. + (-257985. + (116193. - 21126. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (1. + (-980. + (280. + (-63090. + (274015. + (-482475. + (429975. + (-193655. + 35210. * y) * y) * y) * y) * y) * y * y) * y * y) * y * y / 720.),
+ lambda y, s: s * (540. + (540. + (-195. + (-195. + (15. + (63085. + (-274015. + (482475. + (-429975. + (193655. - 35210. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-108. + (-54. + (120. + (60. + (-12. + (-37848. + (164409. + (-289485. + (257985. + (-116193. + 21126. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (12. + (4. + (-15. + (-5. + (3. + (12615. + (-54803. + (96495. + (-85995. + (38731. - 7042. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-1802. + (7829. + (-13785. + (12285. + (-5533. + 1006. * y) * y) * y) * y) * y) * y * y * y * y * y * y / 720.
+ ]
+ elif self.name == 'l6star_c6':
+ self._shift = 3
+ self._weights = [
+ lambda y, s: s * (-12. + (4. + (15. + (-5. + (-3. + (1. + (6587. + (-34869. + (77815. + (-93577. + (63866. + (-23426. + 3604. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (108. + (-54. + (-120. + (60. + (12. + (-6. + (-46109. + (244083. + (-544705. + (655039. + (-447062. + (163982. - 25228. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-540. + (540. + (195. + (-195. + (-15. + (15. + (138327. + (-732249. + (1634115. + (-1965117. + (1341186. + (-491946. + 75684. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (1. + (-980. + (280. + (-20. + (-230545. + (1220415. + (-2723525. + (3275195. + (-2235310. + (819910. - 126140. * y) * y) * y) * y) * y) * y) * y) * y * y) * y * y) * y * y / 720.),
+ lambda y, s: s * (540. + (540. + (-195. + (-195. + (15. + (15. + (230545. + (-1220415. + (2723525. + (-3275195. + (2235310. + (-819910. + 126140. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-108. + (-54. + (120. + (60. + (-12. + (-6. + (-138327. + (732249. + (-1634115. + (1965117. + (-1341186. + (491946. - 75684. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (12. + (4. + (-15. + (-5. + (3. + (1. + (46109. + (-244083. + (544705. + (-655039. + (447062. + (-163982. + 25228. * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y) * y / 720.,
+ lambda y, s: s * (-6587. + (34869. + (-77815. + (93577. + (-63866. + (23426. - 3604. * y) * y) * y) * y) * y) * y) * y * y * y * y * y * y * y / 720.
+ ]
+ else:
+ raise ValueError("Unknown formula : " + self.name)
def slice_i_along_d(self, i, d):
l = list(self._slice_all)
l[d] = i
return tuple(l)
- @abstractmethod
- def __call__(self, ppos, pscal, d):
+ def __call__(self, ppos, origin, pscal, d, res=None,
+ i_row=None, i_col=None, i_dep=None, i_y=None, work=None):
+
"""
Remesh particles at position ppos with scalar
pscal along direction d.
@@ -32,5 +144,52 @@ class Remeshing(NumMethod):
@param ppos : particle position
@param pscal : particle scalar
@param d : splitting direction
+ @param i_row : allocated working array (rows indices array).
+ @param i_col : allocated working array (columns indices array).
+ @param i_dep : allocated working array (depth indices array).
+ @param i_y : allocated working array (distance to left grid point).
"""
- pass
+ shape = ppos.shape
+ if len(shape) == 2:
+ index = tuple([i_row, i_col])
+ if d == 0:
+ index[1][...] = np.indices((shape[1],))[0][np.newaxis, :]
+ else:
+ index[0][...] = np.indices((shape[0],))[0][:, np.newaxis]
+ elif len(shape) == 3:
+ index = tuple([i_row, i_col, i_dep])
+ if d == 0:
+ index[1][...] = \
+ np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
+ index[2][...] = \
+ np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
+ elif d == 1:
+ index[0][...] = \
+ np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
+ index[2][...] = \
+ np.indices((shape[2],))[0][np.newaxis, np.newaxis, :]
+ else:
+ index[0][...] = \
+ np.indices((shape[0],))[0][:, np.newaxis, np.newaxis]
+ index[1][...] = \
+ np.indices((shape[1],))[0][np.newaxis, :, np.newaxis]
+ else:
+ index = tuple([i_row, ])
+
+ floor = res
+ floor[...] = (ppos - origin) / self.dx[d]
+ i_y[...] = (floor)
+ floor[...] = np.floor(floor)
+ i_y[...] -= floor
+
+ index[d][...] = (floor.astype(np.int32) - self._shift) % shape[d]
+ res[...] = 0.
+
+ for w_id, w in enumerate(self._weights):
+ if w_id > 0:
+ index[d][...] = (index[d] + 1) % shape[d]
+ for i in xrange(shape[d]):
+ sl = self.slice_i_along_d(i, d)
+ work = w(i_y, pscal)
+ index_sl = tuple([ind[sl] for ind in index])
+ res[index_sl] = res[index_sl] + work[sl]
diff --git a/HySoP/hysop/operator/discrete/particle_advection.py b/HySoP/hysop/operator/discrete/particle_advection.py
index 44a7a9a57b8d49f08372ac9a653e33ec4c6ac2ee..dc4dc69c928c77f699da39766b0ef0e800b5260b 100644
--- a/HySoP/hysop/operator/discrete/particle_advection.py
+++ b/HySoP/hysop/operator/discrete/particle_advection.py
@@ -28,7 +28,9 @@ class ParticleAdvection(DiscreteOperator):
@param method : Method used
- 'rk2' : Runge Kutta 2nd order advection
- 'm4prime' : M4prime formula
- - 'm6prime' : M6prime formula
+ - 'l4star' : Labmda4star formula
+ - 'l4star_c3' : Labmda4star C3 formula
+ - 'l4star_c4' : Labmda4star C4 formula
- 'm8prime' : M8prime formula
- 'l6star' : Lambda6star formula
"""
@@ -51,6 +53,10 @@ class ParticleAdvection(DiscreteOperator):
self.part_position = part_position
self.part_advectedFields = part_advectedFields
self.num_method = None
+ self.available_methods = [
+ 'm4prime', 'm8prime',
+ 'l4star_c4', 'l4star_c3', 'l4star',
+ 'l6star_c6', 'l6star_c5', 'l6star_c4','l6star']
@debug
def setUp(self):
@@ -75,25 +81,19 @@ class ParticleAdvection(DiscreteOperator):
self.name += '_rk2'
from parmepy.numerics.integrators.runge_kutta2 import RK2_scalar
self.num_advec = RK2_scalar()
- if self.method.find('m4prime') >= 0:
- self.name += '_m4prime'
- from parmepy.numerics.remeshing.m4prime import M4Prime as Remesh
- elif self.method.find('m6prime') >= 0:
- self.name += '_m6prime'
- from parmepy.numerics.remeshing.m6prime import M6Prime as Remesh
- elif self.method.find('m8prime') >= 0:
- self.name += '_m8prime'
- from parmepy.numerics.remeshing.m8prime import M8Prime as Remesh
- elif self.method.find('l6star') >= 0:
- self.name += '_l6star'
- from parmepy.numerics.remeshing.l6star import L6Star as Remesh
- else:
- self.name += '_m6prime'
- from parmepy.numerics.remeshing.m6prime import M6Prime as Remesh
- self.num_remesh = Remesh(
+ remesh = 'l4star'
+ for m in self.available_methods:
+ if self.method.find(m) >= 0:
+ remesh = m
+ break
+ self.name += '_' + remesh
+ from parmepy.numerics.remeshing.remeshing import Remeshing
+ self.num_remesh = Remeshing(
self.velocity.dimension,
self.advectedFields[0].topology.mesh.space_step,
- self.advectedFields[0].topology.mesh.origin)
+ self.advectedFields[0].topology.mesh.origin,
+ remesh)
+
self.interpolation = [
Linear(self.velocity, d,
self.advectedFields[0].topology.mesh.space_step,
@@ -108,9 +108,12 @@ class ParticleAdvection(DiscreteOperator):
self.i_row, self.i_col = \
np.indices(self.part_position.data.shape)
self.i_dep = None
- else:
+ elif len(self.part_position.data.shape) == 3:
self.i_row, self.i_col, self.i_dep = \
np.indices(self.part_position.data.shape)
+ else:
+ self.i_row = np.indices(self.part_position.data.shape)[0]
+ self.i_col, self.i_dep = None, None
self._isUpToDate = True
@@ -152,7 +155,7 @@ class ParticleAdvection(DiscreteOperator):
self.num_advec(
self.interpolation[self.dir],
t, dt, self.part_position.data, res=self.part_position.data,
- work=self.work,
+ work=self.work, yprime=self.velocity.data[self.dir],
i_row=self.i_row, i_col=self.i_col, i_dep=self.i_dep,
i_y=self.i_y)
@@ -166,14 +169,14 @@ class ParticleAdvection(DiscreteOperator):
adF.topology.mesh.origin[self.dir],
p_adF.data[dim], self.dir, res=adF[dim],
i_row=self.i_row, i_col=self.i_col, i_dep=self.i_dep,
- i_y=self.i_y)
+ i_y=self.i_y, work=self.work)
else:
self.num_remesh(
self.part_position.data,
adF.topology.mesh.origin[self.dir],
p_adF.data, self.dir, res=adF.data,
i_row=self.i_row, i_col=self.i_col, i_dep=self.i_dep,
- i_y=self.i_y)
+ i_y=self.i_y, work=self.work)
@debug
def finalize(self):
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index feb952200bd48ba732035b252e7829a629fdcca9..78abbcd43db4056021e3ae98502ad6634458c0a7 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -42,7 +42,7 @@ class Printer(Monitoring):
self.output = []
## Method to collect data in case of distributed data
self.get_data_method = None
- if self.freq != 0:
+ if self.freq > 0:
## Printer core method
self.step = self._printStep
else:
@@ -62,7 +62,14 @@ class Printer(Monitoring):
if simulation is None:
raise ValueError("Missing simulation value for monitoring.")
- self.step(simulation.currentIteration)
+ self._call_number += 1
+ if self.freq > 0 and (simulation.currentIteration % self.freq) == 0:
+ self.step(simulation.currentIteration)
+
+ def finalize(self):
+ if self.freq == -1:
+ self._printStep(self._call_number)
+ Monitoring.finalize(self)
def _build_vtk_dict(self):
"""Build a dictionary from fields to VTK image."""
@@ -105,77 +112,92 @@ class Printer(Monitoring):
If fails, turn to classical ascii output.
@param iter : current iteration number
"""
- if (ite % self.freq) == 0:
- self._call_number += 1
- # Transfer from GPU to CPU if required
- for f in self.variables:
- for df in f.discreteFields.values():
- try:
- df.toHost()
- except AttributeError:
- pass
- # Set output file name
- filename = self.prefix + str(main_rank)
- filename += '_' + "iter_{0:03d}".format(ite) + self.ext
- print "Print to file " + filename
- ## VTK output \todo: Need fix in 2D, getting an IOError.
- if self.ext == '.vtk':
- #orig = self.variables[0].discreteFields.values()[0].topology.mesh.global_start
- topo = self.variables[0].discreteFields.keys()[0]
- orig = tuple([topo.mesh.origin[i]
- for i in xrange(topo.mesh.dim)])
- coords = topo.mesh.coords
- ind = topo.mesh.local_start
- orig = tuple([topo.mesh.coords[i].flatten()[ind[i]]
- for i in xrange(topo.mesh.dim)])
- spacing = tuple(topo.mesh.space_step)
- evtk.imageToVTK(filename, origin=orig, spacing=spacing,
- pointData=self._build_vtk_dict())
- elif self.ext == '.dat':
- if self.variables[0].domain.dimension == 2:
- ## Standard output
- f = open(filename, 'w')
- df = self.variables[0].discreteFields.values()[0]
- shape = df.topology.mesh.resolution
- coords = df.topology.mesh.coords
- if len(shape) == 2:
- for i in xrange(shape[0] - 1):
- for j in xrange(shape[1] - 1):
- f.write("{0:8.12} {1:8.12} ".format(
- coords[0][i, 0], coords[1][0, j]))
- for field in self.variables:
- df = field.discreteFields.values()[0]
- if field.isVector:
- f.write("{0:8.12} {1:8.12} ".format(
- df[0][i, j],
- df[1][i, j]))
- else:
- f.write("{0:8.12} ".format(
- df[i, j]))
- f.write("\n")
- else:
- for i in xrange(shape[0] - 1):
- for j in xrange(shape[1] - 1):
- for k in xrange(shape[2] - 1):
+ # Transfer from GPU to CPU if required
+ for f in self.variables:
+ for df in f.discreteFields.values():
+ try:
+ df.toHost()
+ except AttributeError:
+ pass
+ # Set output file name
+ filename = self.prefix + str(main_rank)
+ filename += '_' + "iter_{0:03d}".format(ite) + self.ext
+ print "Print to file " + filename
+ ## VTK output \todo: Need fix in 2D, getting an IOError.
+ if self.ext == '.vtk':
+ #orig = self.variables[0].discreteFields.values()[0].topology.mesh.global_start
+ topo = self.variables[0].discreteFields.keys()[0]
+ orig = tuple([topo.mesh.origin[i]
+ for i in xrange(topo.mesh.dim)])
+ coords = topo.mesh.coords
+ ind = topo.mesh.local_start
+ orig = tuple([topo.mesh.coords[i].flatten()[ind[i]]
+ for i in xrange(topo.mesh.dim)])
+ spacing = tuple(topo.mesh.space_step)
+ evtk.imageToVTK(filename, origin=orig, spacing=spacing,
+ pointData=self._build_vtk_dict())
+ elif self.ext == '.dat':
+ f = open(filename, 'w')
+ if self.variables[0].domain.dimension == 2:
+ ## Standard output
+ df = self.variables[0].discreteFields.values()[0]
+ shape = df.topology.mesh.resolution
+ coords = df.topology.mesh.coords
+ if len(shape) == 2:
+ for i in xrange(shape[0] - 1):
+ for j in xrange(shape[1] - 1):
+ f.write("{0:8.12} {1:8.12} ".format(
+ coords[0][i, 0], coords[1][0, j]))
+ for field in self.variables:
+ df = field.discreteFields.values()[0]
+ if field.isVector:
+ f.write("{0:8.12} {1:8.12} ".format(
+ df[0][i, j],
+ df[1][i, j]))
+ else:
+ f.write("{0:8.12} ".format(
+ df[i, j]))
+ f.write("\n")
+ elif self.variables[0].domain.dimension == 3:
+ df = self.variables[0].discreteFields.values()[0]
+ shape = df.topology.mesh.resolution
+ coords = df.topology.mesh.coords
+ for i in xrange(shape[0] - 1):
+ for j in xrange(shape[1] - 1):
+ for k in xrange(shape[2] - 1):
+ f.write(
+ "{0:8.12} {1:8.12} {2:8.12} ".format(
+ coords[0][i, 0, 0],
+ coords[1][0, j, 0],
+ coords[2][0, 0, k]))
+ for field in self.variables:
+ df = field.discreteFields.values()[0]
+ if field.isVector:
f.write(
- "{0:8.12} {1:8.12} {2:8.12} ".format(
- coords[0][i, 0, 0],
- coords[1][0, j, 0],
- coords[2][0, 0, k]))
- for field in self.variables:
- df = field.discreteFields.values()[0]
- if field.isVector:
- f.write(
- "{0:8.12} {1:8.12} " +
- "{2:8.12} ".format(
- df[0][i, j, k],
- df[1][i, j, k],
- df[2][i, j, k]))
- else:
- f.write("{0:8.12} ".format(
- df[i, j, k]))
- f.write("\n")
- f.close()
+ "{0:8.12} {1:8.12} " +
+ "{2:8.12} ".format(
+ df[0][i, j, k],
+ df[1][i, j, k],
+ df[2][i, j, k]))
+ else:
+ f.write("{0:8.12} ".format(
+ df[i, j, k]))
+ f.write("\n")
+ else:
+ df = self.variables[0].discreteFields.values()[0]
+ shape = df.topology.mesh.resolution
+ coords = df.topology.mesh.coords
+ for i in xrange(shape[0] - 1):
+ f.write("{0:8.12} ".format(coords[0][i]))
+ for field in self.variables:
+ df = field.discreteFields.values()[0]
+ if field.isVector:
+ f.write(
+ "{0:8.12} ".format(df[0][i]))
+ else:
+ f.write("{0:8.12} ".format(df[i]))
+ f.write("\n")
+ f.close()
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/problem/simulation.py b/HySoP/hysop/problem/simulation.py
index 746180989d50ae7e075f537f4d78cec06ed2720c..b223b9d93b74b23148721f95afd54f2fdb93cb19 100644
--- a/HySoP/hysop/problem/simulation.py
+++ b/HySoP/hysop/problem/simulation.py
@@ -42,6 +42,9 @@ class Simulation(object):
if self.currentIteration < self.iterMax:
self.currentIteration += 1
self.time += self.timeStep
+ # adjust last timestep to reach self.end
+ if self.end - self.time < self.timeStep and self.end > self.time:
+ self.timeStep = self.end - self.time
else:
self.isOver = True
if self.time >= self.end:
@@ -52,8 +55,9 @@ class Simulation(object):
print current state
"""
if main_rank == 0:
- print "=== Iteration : {0:3d}, start from t ={1:6.3f} ===".format(
- self.currentIteration, self.time)
+ print "=== Iteration : {0:3d}, start from t ={1:6.3f} \
+to t ={2:6.3f} ===".format(
+ self.currentIteration, self.time, self.time + self.timeStep)
def reset(self):
"""