diff --git a/docs/sphinx/users_guide/remeshing_and_interpolation.rst b/docs/sphinx/users_guide/remeshing_and_interpolation.rst
index 723586743c63e113d6d4af802488f77a10051a41..a5d69e721a8901b367dded4e11d4d3c5a9dc3398 100644
--- a/docs/sphinx/users_guide/remeshing_and_interpolation.rst
+++ b/docs/sphinx/users_guide/remeshing_and_interpolation.rst
@@ -4,6 +4,14 @@
Interpolation and remeshing methods
===================================
+Interpolation
+-------------
+
+Use :class:`~hysop.numerics.interpolation.Linear to interpolate a field (scalar or vector)
+from one mesh to another.
+
+
+:class:`~hysop.numerics.interpolation.Interpolation
**What**: interpolation of a field, defined as a numpy array from one "grid" to another.
@@ -16,7 +24,11 @@ To interpolate a field, the following inputs are required:
+* requirements for position parameter
+ * :math:`shape[d] == topo_target.mesh.resolution[d] \forall d \neq idir` where idir is the interpolation direction.
+ * :math:`shape[idir] \leq topo_target.mesh.compute_resolution[d]`
+
Devel notes
diff --git a/hysop/mpi/topology.py b/hysop/mpi/topology.py
index ed8dd66856d6dedd9f2e51903ab2dc8b53127ab6..1ec4759fa007d7b91975c7010173db1067183d23 100644
--- a/hysop/mpi/topology.py
+++ b/hysop/mpi/topology.py
@@ -155,9 +155,17 @@ class Cartesian(object):
# previous (resp. next) neighbour in direction i
# (warning : direction in the grid of process).
self.neighbours = npw.dim_zeros((2, self.dimension))
+ self.neighbours_tmp = []
+ for direction in range(self.domain.dimension):
+ self.neighbours_tmp.append(self.comm.Shift(direction, 1))
+ self.distributed_dirs = tuple([d for d in xrange(self.domain.dimension)
+ if self.neighbours_tmp[d][0]
+ is not self.rank])
+ self.bc_dirs = tuple([d for d in xrange(self.domain.dimension)
+ if d not in self.distributed_dirs])
+
for direction in range(self.dimension):
self.neighbours[:, direction] = self.comm.Shift(direction, 1)
-
# ===== 4 - Computation of the local meshes =====
# Local mesh on the current mpi process.
@@ -475,6 +483,8 @@ class Cartesian(object):
s += str(self._mpis.rank) + '.\n'
s += '- neighbours coordinates : \n'
s += str(self.neighbours) + '\n'
+ s += '- FULL neighbours coordinates : \n'
+ s += str(self.neighbours_tmp) + '\n'
s += str(self.mesh)
s += '\n=================================\n'
return s
diff --git a/hysop/numerics/interpolation.py b/hysop/numerics/interpolation.py
index b77f333da4ede19c20b504027158650ab11a9ddb..505d0a4ae0a097394e587357413f9d2b9284c17e 100644
--- a/hysop/numerics/interpolation.py
+++ b/hysop/numerics/interpolation.py
@@ -190,9 +190,9 @@ class Interpolation(object):
except in direction of interpolation where any number of points
is allowed, with a max of points equal to the number of grid
points in this direction (excluding ghosts).
- work : list of numpy arrays.
- Used as local buffer. work[0] must be of the same shape
- as positions. work[0] contains iy in return,
+ work : numpy array.
+ Used as local buffer. work must be of the same shape
+ as positions and contains iy in return,
i_y = (xp - xg)/dx
Returns
@@ -212,15 +212,16 @@ class Interpolation(object):
"""
# local origin and space step in current direction
x0 = self.topo_source.mesh.origin[self.direction]
- i_y = work[0]
- assert npw.arrays_share_data(i_y, self._rwork[0])
+ #i_y = work
+ # memory optim: check if input work 'belongs' to rwork
+ #assert npw.arrays_share_data(i_y, work)
# -- set indices of closest grid points, on the left --
# Initialized with source grid coordinates
left_points_indices = list(self._icoords)
# = [self._icoords[i] for i in xrange(dimension)]
# pick memory from self._iwork
iwork_dir = WorkSpaceTools.check_work_array(
- 1, i_y[...].shape, self._iwork, data_type=HYSOP_INTEGER)
+ 1, work.shape, self._iwork, data_type=HYSOP_INTEGER)
# we need an array (filled later) for current direction
left_points_indices[self.direction] = iwork_dir[0]
# check memory consumption
@@ -228,13 +229,13 @@ class Interpolation(object):
self._iwork[0])
# -- compute distance between particles positions and
# closest grid point, on the left --
- i_y[...] = (positions[self._target_indices] - x0) *\
+ work[...] = (positions[self._target_indices] - x0) *\
self._inv_dx[self.direction]
#self._assert_particles_roundoff(work, self._target_indices)
# update interpolation point indices (in self._iwork)
left_points_indices[self.direction][...] = \
- npw.asintarray(np.floor(i_y[...]))
- i_y[...] -= left_points_indices[self.direction][...]
+ npw.asintarray(np.floor(work[...]))
+ work[...] -= left_points_indices[self.direction][...]
nmax = self.topo_source.mesh.resolution[self.direction]
# Ensure all points belongs to local domain
# (including ghosts)
@@ -474,8 +475,8 @@ class Linear(Interpolation):
Notes
-----
* since interpolation function is likely to be used
- as right-hand side of an ode solver, the argument slist must fit with
- the integrator requirements, which explains the 't' argument
+ as right-hand side of an ode solver, the list of arguments must fit
+ with the integrator requirements, which explains the 't' argument
unused below but required, or the positions and result as lists.
"""
assert isinstance(positions, list)
@@ -500,16 +501,24 @@ class Linear(Interpolation):
# pick buffer memory from self._rwork
# wk[0] shares memory with rwork[0] but data are rearranged
# according to positions[0].shape.
+ print nb_points, self.direction
wk = WorkSpaceTools.check_work_array(
1, positions[0][self._target_indices].shape, self._rwork)
- wk.append(result[0])
- left_points_indices = self._compute_iy(positions[0], wk)
+ # memory optim: check if input work 'belongs' to rwork
+ assert npw.arrays_share_data(self._rwork[0], wk[0])
+ left_points_indices = self._compute_iy(positions[0], wk[0])
assert npw.arrays_share_data(left_points_indices[self.direction],
self._iwork[0])
- wk[1][self._target_indices] = \
+ print self.interpolated_field[left_points_indices].shape
+ #refsize = positions[0].size
+ #refshape = positions[0].shape
+ #wk.append(result[0][])
+ #res = WorkSpaceTools.check_work_array(1, refshape, result)
+ result[0][self._target_indices] = \
self.interpolated_field[left_points_indices] * (1. - wk[0][...])
# compute 'right points' indices
left_points_indices[self.direction] += 1
- wk[1][self._target_indices] += \
+
+ result[0][self._target_indices] += \
self.interpolated_field[left_points_indices] * wk[0][...]
return result
diff --git a/hysop/numerics/odesolvers.py b/hysop/numerics/odesolvers.py
index 33ab3027f8ba17c23dd58698277fa8252c1b19b0..79741682d402582df901f69bc959aafebcfe3547 100755
--- a/hysop/numerics/odesolvers.py
+++ b/hysop/numerics/odesolvers.py
@@ -31,7 +31,7 @@ class ODESolver(object):
__metaclass__ = ABCMeta
def __init__(self, nb_components, topo, f=None, indices=None,
- rwork=None, optim=None):
+ rwork=None, optim=None, wk_shape=None, need_synchro=False):
"""
Parameters
----------
@@ -56,6 +56,13 @@ class ODESolver(object):
optim : int, optional
the level of optimization (memory management).
Default = None.
+ wk_shape: tuple, optional
+ shape of the internal work vector. Default = input topo
+ resolution.
+ need_synchro: boolean, optional
+ Set to true if y parameter of the rhs needs an update
+ of its ghost points before each call. This might be
+ the case if rhs is an interpolation. Default=False
Notes
-----
@@ -65,6 +72,9 @@ class ODESolver(object):
system is solved. This is the default behavior
* optim = WITH_GUESS : in that case, work must contain a first
evaluation of f(t,y) before each call.
+ * wk_shape might be usefull when the advected field shape
+ in advection dir is lower than topo shape (for instance when
+ pushing particles with a threshhold parameter).
"""
# RHS.
@@ -88,19 +98,21 @@ class ODESolver(object):
#self.output_indices = None
# work arrays
- self.rwork = self._set_work_arrays(topo, rwork=rwork)
- self._synchronize = UpdateGhosts(topo, self._nb_components)
+ self.rwork = self._set_work_arrays(topo, wk_shape, rwork=rwork)
+ if need_synchro:
+ self._synchronize = UpdateGhosts(topo, self._nb_components)
+ self.compute_rhs = self._update_rhs_with_synchro
+ else:
+ self.compute_rhs = self.f
- def _set_work_arrays(self, topo, reduce_output_shape=False, rwork=None):
+ def _set_work_arrays(self, topo, wk_shape=False, rwork=None):
"""Check and allocate internal work buffers.
"""
wk_prop = self.get_work_properties(self._nb_components, topo)['rwork']
if wk_prop is None:
return []
- if reduce_output_shape:
- subsize = np.asarray([self.in_indices[i].stop -
- self.in_indices[i].start
- for i in xrange(len(self.in_indices))])
+ if wk_shape:
+ subsize = wk_shape
else:
subsize = topo.mesh.resolution
@@ -118,6 +130,19 @@ class ODESolver(object):
topology on which the integrator will be applied.
"""
+ # def update_rhs_local(self, t, y, result):
+ # """Return rhs value, no synchronization between mpi domains.
+ # """
+ # return self.f(t, y, result)
+
+ def _update_rhs_with_synchro(self, t, y, result):
+ """"Return rhs value, but first
+ synchronize y values between mpi domains
+ (i.e. ghost points update)
+ """
+ self._synchronize(y)
+ return self.f(t, y, result)
+
def __call__(self, t, y, dt, result):
"""Apply integrator
@@ -161,7 +186,7 @@ class ODESolver(object):
result may be equal to y.
"""
self.rwork[:self._nb_components] = \
- self.f(t, y, self.rwork[:self._nb_components])
+ self.compute_rhs(t, y, self.rwork[:self._nb_components])
return self._core(t, y, dt, result)
@abstractmethod
@@ -295,7 +320,7 @@ class RK2(ODESolver):
Note : since result may be equal to y, it can not
be used as a temporary workspace.
"""
- ic = self.in_indices
+ i_y = self.in_indices
for i in xrange(self._nb_components):
for j in xrange(len(self.rwork)):
assert not npw.arrays_share_data(y[i], self.rwork[j])
@@ -308,19 +333,20 @@ class RK2(ODESolver):
# result = y + dt * k2
# yn
+ #print ic, y[0].shape
for i in xrange(self._nb_components):
- np.add(y[i][ic], work0[i] * 0.5 * dt, yn[i])
+ np.add(y[i][i_y], work0[i] * 0.5 * dt, yn[i])
+ #np.add(y[i][ic], work0[i] * 0.5 * dt, yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# k2 in work0
- work0 = self.f(t + 0.5 * dt, yn, work0)
+ work0 = self.compute_rhs(t + 0.5 * dt, yn, work0)
# *= dt
for i in xrange(self._nb_components):
#np.multiply(work0[i], dt, work0[i])
work0[i] *= dt
# result = y + work0
- np.add(work0[i], y[i][ic], result[i])
+ np.add(work0[i], y[i][i_y], result[i][i_y])
+ #np.add(work0[i], y[i][ic], result[i])
return result
@@ -396,18 +422,14 @@ class RK3(ODESolver):
for i in xrange(self._nb_components):
np.add(y[i], work0[i] * dt / 3, yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# k2 in kn
- kn = self.f(t + dt / 3, yn, kn)
+ kn = self.compute_rhs(t + dt / 3, yn, kn)
# yn
for i in xrange(self._nb_components):
np.add(y[i], 2 * dt / 3 * kn[i], yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# k3 in kn
- kn = self.f(t + 2 * dt / 3, yn, kn)
+ kn = self.compute_rhs(t + 2 * dt / 3, yn, kn)
# k1 + 3 * k3 in work0
for i in xrange(self._nb_components):
np.add(work0[i], 3 * kn[i], work0[i])
@@ -483,10 +505,8 @@ class RK4(ODESolver):
for i in xrange(self._nb_components):
np.add(y[i], work0[i] * dt / 2, yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# k2 in kn
- kn = self.f(t + dt / 2, yn, kn)
+ kn = self.compute_rhs(t + dt / 2, yn, kn)
# k1 + 2 * k2 in work0
for i in xrange(self._nb_components):
@@ -494,10 +514,8 @@ class RK4(ODESolver):
# yn
np.add(y[i], dt / 2 * kn[i], yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# k3 in kn
- kn = self.f(t + dt / 2, yn, kn)
+ kn = self.compute_rhs(t + dt / 2, yn, kn)
# k1 + 2 * k2 + 2 * k3 in work0
for i in xrange(self._nb_components):
@@ -505,10 +523,8 @@ class RK4(ODESolver):
# yn
np.add(y[i], dt * kn[i], yn[i])
- # Update ghosts
- self._synchronize.apply(yn)
# K4 in kn
- kn = self.f(t + dt, yn, kn)
+ kn = self.compute_rhs(t + dt, yn, kn)
# k1 + 2 * k2 + 2 * k3 + k4
for i in xrange(self._nb_components):
diff --git a/hysop/numerics/remeshing.py b/hysop/numerics/remeshing.py
index f63f52ecad8ad3c12fd7d33f80d78f27e0430a61..80e17ece4ffcf9840b7c85a27da4bacbfd9f1aa5 100644
--- a/hysop/numerics/remeshing.py
+++ b/hysop/numerics/remeshing.py
@@ -114,7 +114,8 @@ class Remeshing(Interpolation):
# wk[0] --> iy, wk[1] --> used in weights accumulation below
wk = WorkSpaceTools.check_work_array(
2, ppos[self._target_indices].shape, self._rwork)
- left_points_indices = self._compute_iy(ppos, wk)
+ assert npw.arrays_share_data(self._rwork[0], wk[0])
+ left_points_indices = self._compute_iy(ppos, wk[0])
assert npw.arrays_share_data(left_points_indices[self.direction],
self._iwork[0])
# shift iwork, so that iwork[direction][p] is the index
diff --git a/hysop/numerics/tests/test_update_ghosts.py b/hysop/numerics/tests/test_update_ghosts.py
index b4cfaf2ee219b5729c2d26d594a25338f51f99a2..1a4a7fa9be3d7ccb097836593749d25fb7c03354 100755
--- a/hysop/numerics/tests/test_update_ghosts.py
+++ b/hysop/numerics/tests/test_update_ghosts.py
@@ -1,154 +1,267 @@
+"""Test ghost points synchronization
+"""
from hysop.numerics.update_ghosts import UpdateGhosts, UpdateGhostsFull
from hysop.domain.box import Box
from hysop.fields.continuous import Field
from hysop.tools.parameters import Discretization
+from hysop.mpi import main_rank
+import hysop.tools.numpywrappers as npw
import numpy as np
-def _setup(ghosts, topo_dim):
- dom = Box()
-
+def get_gh_slices(topo):
+ """Return indices of ghost points"""
+ tab = Field(topo.domain, name='temp')
+ tabd = tab.discretize(topo)[0]
+ tabd[...] = 0.
+ tabd[topo.mesh.compute_index] = 1.
+ return np.where(tabd == 0.)
+
+
+def get_ghosts_indices(topo):
+ """Return indices of points inside the ghost layers
+ for a given topology
+
+ gh_s[d] : for points at the beginning of the domain
+ in direction d.
+ gh_e[d] : for points at the end of the domain
+ in direction d.
+
+ Corner points are not included.
+ """
+ gh_s = []
+ gh_e = []
+ ghosts = topo.ghosts()
+ sl_start = [slice(0, ghosts[d]) for d in xrange(topo.domain.dimension)]
+ sl_end = [slice(-ghosts[d], None, None)
+ for d in xrange(topo.domain.dimension)]
+ for j in xrange(topo.domain.dimension):
+ gh_s.append(list(topo.mesh.compute_index))
+ gh_e.append(list(topo.mesh.compute_index))
+ gh_s[j][j] = sl_start[j]
+ gh_e[j][j] = sl_end[j]
+ return gh_s, gh_e
+
+
+def get_corners_indices(topo):
+ """Return indices of points inside the ghost layers
+ for a given topology
+
+ gh_s[d] : for points at the beginning of the domain
+ in direction d.
+ gh_e[d] : for points at the end of the domain
+ in direction d.
+
+ Corner points are included.
+ """
+ gh_s = []
+ gh_e = []
+ ghosts = topo.ghosts()
+ sl_start = [slice(0, ghosts[d]) for d in xrange(topo.domain.dimension)]
+ sl_end = [slice(-ghosts[d], None, None)
+ for d in xrange(topo.domain.dimension)]
+ for j in xrange(topo.domain.dimension):
+ gh_s.append(list(topo.mesh.local_index))
+ gh_e.append(list(topo.mesh.local_index))
+ gh_s[j][j] = sl_start[j]
+ gh_e[j][j] = sl_end[j]
+
+ start, end = get_ghosts_indices(topo)
+ tab = np.zeros((topo.mesh.resolution), dtype=np.int16)
+ tab[topo.mesh.compute_index] = 1
+ for d in xrange(topo.domain.dimension):
+ tab[start[d]] = 1
+ tab[end[d]] = 1
+ corners = np.where(tab == 0)
+
+ corners = [list(topo.mesh.local_index), ] * (2 ** (topo.domain.dimension))
+
+ dimension = topo.domain.dimension
+ corners_size = 2 ** dimension
+ corners = [None, ] * corners_size
+ val = npw.int_zeros(corners_size)
+ if dimension == 1:
+ return [], 0.
+ elif dimension == 2:
+ corners[0] = [sl_start[0], sl_start[1]]
+ corners[1] = [sl_start[0], sl_end[1]]
+ corners[2] = [sl_end[0], sl_start[1]]
+ corners[3] = [sl_end[0], sl_end[1]]
+ val[0] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1, topo.proc_coords[1] - 1])
+ val[1] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1, topo.proc_coords[1] + 1])
+ val[2] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1, topo.proc_coords[1] - 1])
+ val[3] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1, topo.proc_coords[1] + 1])
+ elif dimension == 3:
+ corners[0] = [sl_start[0], sl_start[1], sl_start[2]]
+ corners[1] = [sl_start[0], sl_start[1], sl_end[2]]
+ corners[2] = [sl_start[0], sl_end[1], sl_start[2]]
+ corners[3] = [sl_start[0], sl_end[1], sl_end[2]]
+ corners[4] = [sl_end[0], sl_start[1], sl_start[2]]
+ corners[5] = [sl_end[0], sl_start[1], sl_end[2]]
+ corners[6] = [sl_end[0], sl_end[1], sl_start[2]]
+ corners[7] = [sl_end[0], sl_end[1], sl_end[2]]
+ val[0] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1,
+ topo.proc_coords[1] - 1,
+ topo.proc_coords[2] - 1])
+ val[1] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1,
+ topo.proc_coords[1] - 1,
+ topo.proc_coords[2] + 1])
+ val[2] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1,
+ topo.proc_coords[1] + 1,
+ topo.proc_coords[2] - 1])
+ val[3] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] - 1,
+ topo.proc_coords[1] + 1,
+ topo.proc_coords[2] + 1])
+ val[4] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1,
+ topo.proc_coords[1] - 1,
+ topo.proc_coords[2] - 1])
+ val[5] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1,
+ topo.proc_coords[1] - 1,
+ topo.proc_coords[2] + 1])
+ val[6] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1,
+ topo.proc_coords[1] + 1,
+ topo.proc_coords[2] - 1])
+ val[7] = topo.comm.Get_cart_rank(
+ [topo.proc_coords[0] + 1,
+ topo.proc_coords[1] + 1,
+ topo.proc_coords[2] + 1])
+
+# print "youou ", corners
+ return corners, val
+
+
+def run_synchro(ghosts, cutdir=None, full=False):
+ """Build domain, fields and topology
+ for a given ghost layer size and
+ distributed directions.
+ Then call synchronization and check.
+ """
+ dimension = len(ghosts)
+ dom = Box(length=[1., ] * dimension)
f = Field(dom, is_vector=True, name='f')
- resolTopo = Discretization([33, 33, 33], ghosts=ghosts)
- topo = dom.create_topology(resolTopo, dim=topo_dim)
- f.discretize(topo)
- df = f.discreteFields[topo]
- for i, dfd in enumerate(df.data):
- dfd[...] = 0.
- dfd[topo.mesh.compute_index] = 1. * (i + 1)
- return df, topo
-
-
-def verify(df, gh):
- sli = (
- (slice(0, gh[0]), slice(gh[1], -gh[1]), slice(gh[2], -gh[2])),
- (slice(-gh[0], None), slice(gh[1], -gh[1]), slice(gh[2], -gh[2])),
- (slice(gh[0], -gh[0]), slice(0, gh[1]), slice(gh[2], -gh[2])),
- (slice(gh[0], -gh[0]), slice(-gh[1], None), slice(gh[2], -gh[2])),
- (slice(gh[0], -gh[0]), slice(gh[1], -gh[1]), slice(0, gh[2])),
- (slice(gh[0], -gh[0]), slice(gh[1], -gh[1]), slice(-gh[2], None)))
- slni = (
- (slice(0, gh[0]), slice(0, gh[1]), slice(None)),
- (slice(0, gh[0]), slice(None), slice(0, gh[2])),
- (slice(None), slice(0, gh[1]), slice(0, gh[2])),
- (slice(-gh[0], None), slice(-gh[1], None), slice(None)),
- (slice(-gh[0], None), slice(None), slice(-gh[2], None)),
- (slice(None), slice(-gh[1], None), slice(-gh[2], None)))
- for i, dfd in enumerate(df.data):
- for s in sli:
- assert np.max(dfd[s]) == np.min(dfd[s]) and \
- np.max(dfd[s]) == 1. * (i + 1)
- for s in slni:
- assert np.max(dfd[s]) == np.min(dfd[s]) and \
- np.max(dfd[s]) == 0.
-
-
-def verify_full(df, gh):
- for i, dfd in enumerate(df.data):
- for s in ((slice(0, gh[0]), slice(None), slice(None)),
- (slice(-gh[0], None), slice(None), slice(None)),
- (slice(None), slice(0, gh[1]), slice(None)),
- (slice(None), slice(-gh[1], None), slice(None)),
- (slice(None), slice(None), slice(0, gh[2])),
- (slice(None), slice(None), slice(-gh[2], None))):
- assert np.max(dfd[s]) == np.min(dfd[s]) and \
- np.max(dfd[s]) == 1. * (i + 1)
- assert np.max(dfd) == np.min(dfd) and \
- np.max(dfd) == 1. * (i + 1)
-
-
-def test_update_ghosts_simple_1D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 1)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_simple_2D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 2)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_simple_3D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 3)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_1D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 1)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_2D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 2)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_3D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 3)
- update = UpdateGhosts(topo, 3)
- update(df.data)
- verify(df, gh)
-
-
-def test_update_ghosts_full_simple_1D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 1)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
-
-
-def test_update_ghosts_full_simple_2D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 2)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
-
-
-def test_update_ghosts_full_simple_3D():
- gh = [1, 1, 1]
- df, topo, = _setup(gh, 3)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
-
-
-def test_update_ghosts_full_1D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 1)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
-
-
-def test_update_ghosts_full_2D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 2)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
-
-
-def test_update_ghosts_full_3D():
- gh = [2, 3, 4]
- df, topo, = _setup(gh, 3)
- update = UpdateGhostsFull(topo, 3)
- update(df.data)
- verify_full(df, gh)
+ discr = Discretization([33, ] * dimension, ghosts=ghosts)
+ topo = dom.create_topology(discr, cutdir=cutdir)
+ df = f.discretize(topo)
+ # initialize + 0 in ghost layer
+ nbc = f.nb_components
+ for i in xrange(nbc):
+ df[i][...] = main_rank
+ #df[i][topo.mesh.compute_index] = topo.rank
+
+ for i in xrange(df.nb_components):
+ assert np.allclose(df[i], main_rank)
+
+ if full:
+ sync = UpdateGhostsFull(topo, nbc)
+
+ else:
+ sync = UpdateGhosts(topo, nbc)
+
+ start, end = get_ghosts_indices(topo)
+ corners, val = get_corners_indices(topo)
+
+ sync.apply(df.data)
+ nghb = topo.neighbours_tmp
+ for i in xrange(df.nb_components):
+ assert np.allclose(df[i][topo.mesh.compute_index], main_rank)
+ for d in xrange(topo.domain.dimension):
+ assert np.allclose(df[i][start[d]], nghb[d][0])
+ assert np.allclose(df[i][end[d]], nghb[d][1])
+ if full:
+ for d in xrange(len(corners)):
+ assert np.allclose(df[i][corners[d]], val[d])
+
+
+def test_1_1_simple():
+ """ 1D domain, 1D topo, simple update
+ """
+ run_synchro([1, ])
+
+
+def test_2_1_simple():
+ """ 2D domain, 1D topo, simple update
+ """
+ run_synchro([1, 1], [True, False])
+
+
+def test_2_2_simple():
+ """ 2D domain, 2D topo, simple update
+ """
+ run_synchro([1, 1])
+
+
+def test_3_1_simple():
+ """ 3D domain, 1D topo, simple update
+ """
+ run_synchro([1, 1, 1], [True, False, False])
+
+
+def test_3_2_simple():
+ """ 3D domain, 1D topo, simple update
+ """
+ run_synchro([1, 1, 1], [True, False, True])
+
+
+def test_3_3_simple():
+ """ 3D domain, 1D topo, simple update
+ """
+ run_synchro([1, 1, 1])
+
+
+def test_3_3_bis():
+ """ 3D domain, 3D topo, different ghost size in each dir
+ """
+ run_synchro([1, 3, 2])
+
+
+def test_1_1_full():
+ """ 1D domain, 1D topo, corners update
+ """
+ run_synchro([1, ], full=True)
+
+
+def test_2_1_full():
+ """ 2D domain, 1D topo, corners update
+ """
+ run_synchro([1, 1], [True, False], full=True)
+
+
+def test_2_2_full():
+ """ 2D domain, 2D topo, corners update
+ """
+ run_synchro([1, 1], full=True)
+
+
+def test_3_1_full():
+ """ 3D domain, 1D topo, corners update
+ """
+ run_synchro([1, 1, 1], [True, False, False], full=True)
+
+
+def test_3_2_full():
+ """ 3D domain, 2D topo, corners update
+ """
+ run_synchro([1, 1, 1], [True, False, True], full=True)
+
+
+def test_3_3_full():
+ """ 3D domain, 3D topo, corners update
+ """
+ run_synchro([1, 1, 1], full=True)
+
+def test_3_2_bis_full():
+ """ 3D domain, 2D topo, corners update, different ghost size in each dir
+ """
+ run_synchro([1, 3, 2], [True, False, True], full=True)
diff --git a/hysop/numerics/update_ghosts.py b/hysop/numerics/update_ghosts.py
index 7e6df40c2baca2ad51dc5c7da0d24a994ce42148..4f4a254df92b2bf722536917e7dcaccf77f6a85c 100644
--- a/hysop/numerics/update_ghosts.py
+++ b/hysop/numerics/update_ghosts.py
@@ -60,10 +60,12 @@ class UpdateGhosts(object):
self._g_tonext = []
# Indices of points to be sent to previous neighbour.
self._g_toprevious = []
- # Buffers that contains ghosts points to be sent (list, one per dir)
- self._sendbuffer = []
- # Buffers for reception (list, one per dir)
- self._recvbuffer = []
+ # Buffers that contains ghosts points to be sent
+ self._sendbuffer = None
+ # List of sizes of buffer to be sent/recv in each direction
+ self._buffer_size = []
+ # Buffers for reception
+ self._recvbuffer = None
# domain dimension
self._dim = self.topology.domain.dimension
@@ -75,30 +77,47 @@ class UpdateGhosts(object):
self._memoryblocksize = np.prod(self.memshape)
# Number of numpy arrays that will be updated
self.nb_elements = nb_elements
+
+ # list of functions, _apply[d] corresponds to the synchronization
+ # used in direction d, either
+ # - local (non-distributed) --> apply BC
+ # - distributed --> mpi comm
+ self._apply = [None, ] * self._dim
+ # distributed (mpi) directions
+ self.exchange_dir = []
+ if self.topology.size > 1:
+ self.exchange_dir = [d for d in xrange(self._dim)
+ if self.topology.cutdir[d]]
+ for d in self.exchange_dir:
+ self._apply[d] = self._apply_in_dir
+ for d in self.topology.bc_dirs:
+ self._apply[d] = self._apply_bc_in_dir
+
if self.topology.size > 1: # else no need to set buffers ...
# Size computation below assumes that what we send in one
# dir has the same size as what we receive from process in the
# same dir ...
- exchange_dir = [d for d in xrange(self._dim)
- if self.topology.cutdir[d]]
# A temporary array used to calculate slices sizes
temp = np.zeros(self.memshape, dtype=np.int8)
- for d in exchange_dir:
+ for d in self.exchange_dir:
buffsize = 0
buffsize += temp[self._g_tonext[d]].size * self.nb_elements
- self._sendbuffer.append(npw.zeros((buffsize)))
- self._recvbuffer.append(npw.zeros((buffsize)))
+ self._buffer_size.append(buffsize)
+ self._recvbuffer = npw.zeros(max(self._buffer_size))
+ self._sendbuffer = npw.zeros(max(self._buffer_size))
def _setup_slices(self):
"""Compute slices used to describe what to send
and receive in ghosts points.
-
"""
+ # all points
defslice = [slice(None, None, None)] * self._dim
+ # no points
nogh_slice = [slice(0)] * self._dim
for d in xrange(self._dim):
if self.ghosts[d] > 0:
+ # get ghost points indices in current direction ...
self._g_fromprevious.append(list(defslice))
self._g_fromprevious[d][d] = slice(self.ghosts[d])
self._g_fromnext.append(list(defslice))
@@ -111,6 +130,8 @@ class UpdateGhosts(object):
-self.ghosts[d])
other_dim = self._other_dim(d)
+ # ... and only compute points indices in other directions
+ # i.e. 'corner' points are not taken into account.
for d2 in other_dim:
self._g_fromprevious[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
@@ -121,12 +142,29 @@ class UpdateGhosts(object):
self._g_tonext[d][d2] = slice(self.ghosts[d2],
-self.ghosts[d2])
else:
+ # empty lists of indices
self._g_fromprevious.append(list(nogh_slice))
self._g_fromnext.append(list(nogh_slice))
self._g_toprevious.append(list(nogh_slice))
self._g_tonext.append(list(nogh_slice))
+ self._g_fromprevious = self._immutable(self._g_fromprevious)
+ self._immutable(self._g_fromnext)
+ self._immutable(self._g_toprevious)
+ self._immutable(self._g_tonext)
+
+ @staticmethod
+ def _immutable(var):
+ """Convert list of lists to immutable tuples
+ """
+ assert isinstance(var, list)
+ for i in xrange(len(var)):
+ var[i] = tuple(var[i])
+ return tuple(var)
+
def _other_dim(self, d):
+ """utility to get all directions except d
+ """
return [x for x in xrange(self._dim) if x != d]
def __call__(self, variables):
@@ -145,43 +183,71 @@ class UpdateGhosts(object):
assert (self.topology.domain.boundaries == PERIODIC).all(),\
'Only implemented for periodic boundary conditions.'
assert isinstance(variables, list)
- dirs = [d for d in xrange(self._dim)
- if self.topology.shape[d] == 1]
+ dirs = self.topology.bc_dirs
for d in dirs:
- self._apply_bc_in_dir(variables, d)
+ self._apply_bc_in_dir(variables, d, 0)
- def _apply_bc_in_dir(self, variables, d):
+ def _apply_bc_in_dir(self, variables, d, i):
"""Apply periodic boundary condition in direction d."""
+ assert (self.topology.domain.boundaries == PERIODIC).all(),\
+ 'Only implemented for periodic boundary conditions.'
for v in variables:
assert v.shape == self.memshape
v[self._g_fromprevious[d]] = v[self._g_tonext[d]]
v[self._g_fromnext[d]] = v[self._g_toprevious[d]]
+ return i
+ # @debug
+ # def apply(self, variables):
+ # """Update ghosts points values:
+ # mpi communications and boundary conditions.
+
+ # Parameters
+ # ----------
+ # variables : list of numpy arrays
+ # arrays that must be synchronized.
+ # """
+ # assert isinstance(variables, list)
+ # # Apply boundary conditions for distributed directions ...
+ # i = 0
+ # for d in self.exchange_dir:
+ # self._apply_in_dir(variables, d, i)
+ # # update position in buffers list
+ # i += 1
+ # # ... and for non-distributed directions.
+ # self.apply_bc(variables)
@debug
def apply(self, variables):
- """
- Compute ghosts values from mpi communications and boundary conditions.
+ """Update ghosts points values:
+ mpi communications and boundary conditions.
+
+ Parameters
+ ----------
+ variables : list of numpy arrays
+ arrays that must be synchronized.
"""
assert isinstance(variables, list)
- exchange_dir = []
- if self.topology.size > 1:
- exchange_dir = [d for d in xrange(self._dim)
- if self.topology.cutdir[d]]
+ # Apply boundary conditions all directions ...
i = 0
- for d in exchange_dir:
- self._apply_in_dir(variables, d, i)
- # update index in neighbours list
- i += 1
- # End of loop through send/recv directions.
- # Apply boundary conditions for non-distributed directions
- self.apply_bc(variables)
+ for d in xrange(self._dim):
+ i = self._apply[d](variables, d, i)
def _apply_in_dir(self, variables, d, i):
"""Communicate ghosts values in direction d for neighbour
- in direction i of the topology"""
+ in direction i of the topology
+
+ Parameters
+ ----------
+ variables : list of numpy arrays
+ arrays to be synchronized
+ d : int
+ direction of synchronization
+ i : int
+ counter, current position in buffers.
+ """
comm = self.topology.comm
rank = self.topology.rank
- neighbours = self.topology.neighbours
+ neighbours = self.topology.neighbours_tmp[d]
# 1 - Fill in buffers
# Loop through all variables that are distributed
pos = 0
@@ -189,15 +255,18 @@ class UpdateGhosts(object):
for v in variables:
assert v.shape == self.memshape
nextpos += v[self._g_tonext[d]].size
- self._sendbuffer[i][pos:nextpos] = v[self._g_tonext[d]].flat
+ self._sendbuffer[pos:nextpos] = v[self._g_tonext[d]].flat
pos = nextpos
# 2 - Send to next receive from previous
- dest_rk = neighbours[1, i]
- from_rk = neighbours[0, i]
- comm.Sendrecv([self._sendbuffer[i], HYSOP_MPI_REAL],
+ dest_rk = neighbours[1]
+ from_rk = neighbours[0]
+ sendbuffer = self._sendbuffer[:self._buffer_size[i]]
+ recvbuffer = self._recvbuffer[:self._buffer_size[i]]
+
+ comm.Sendrecv([sendbuffer, HYSOP_MPI_REAL],
dest=dest_rk, sendtag=rank,
- recvbuf=self._recvbuffer[i],
+ recvbuf=recvbuffer,
source=from_rk, recvtag=from_rk)
# 3 - fill variables with recvbuffer and update sendbuffer
@@ -207,17 +276,17 @@ class UpdateGhosts(object):
for v in variables:
nextpos += v[self._g_fromprevious[d]].size
v[self._g_fromprevious[d]].flat = \
- self._recvbuffer[i][pos:nextpos]
- self._sendbuffer[i][pos:nextpos] = \
+ self._recvbuffer[pos:nextpos]
+ self._sendbuffer[pos:nextpos] = \
v[self._g_toprevious[d]].flat
pos = nextpos
# 4 -Send to previous and receive from next
- dest_rk = neighbours[0, i]
- from_rk = neighbours[1, i]
- comm.Sendrecv([self._sendbuffer[i], HYSOP_MPI_REAL],
+ dest_rk = neighbours[0]
+ from_rk = neighbours[1]
+ comm.Sendrecv([sendbuffer, HYSOP_MPI_REAL],
dest=dest_rk, sendtag=rank,
- recvbuf=self._recvbuffer[i],
+ recvbuf=recvbuffer,
source=from_rk, recvtag=from_rk)
# 5 - fill variables with recvbuffer
pos = 0
@@ -225,9 +294,12 @@ class UpdateGhosts(object):
for v in variables:
nextpos += v[self._g_fromprevious[d]].size
v[self._g_fromnext[d]].flat = \
- self._recvbuffer[i][pos:nextpos]
+ self._recvbuffer[pos:nextpos]
pos = nextpos
+ # update position in buffers list
+ return i + 1
+
class UpdateGhostsFull(UpdateGhosts):
"""Ghost points synchronization for a list of numpy arrays
@@ -250,27 +322,22 @@ class UpdateGhostsFull(UpdateGhosts):
# if d == 0:
return [x for x in xrange(self._dim) if x > d]
- @debug
- def apply(self, variables):
- """Apply either mpi communications
- or local boundary conditions to fill ghosts.
- Loop over directions and switch among local BC or mpi comm.
- """
- assert isinstance(variables, list)
- exchange_dir = []
- if self.topology.size > 1:
- exchange_dir = [d for d in xrange(self._dim)
- if self.topology.cutdir[d]]
- local_bc_dir = [d for d in xrange(self._dim)
- if self.topology.shape[d] == 1]
- assert len(exchange_dir) + len(local_bc_dir) == self._dim
-
- i = 0
- for d in xrange(self._dim):
- if d in local_bc_dir:
- self._apply_bc_in_dir(variables, d)
- elif d in exchange_dir:
- self._apply_in_dir(variables, d, i)
- # update index in neighbours list
- i += 1
- # End of loop through send/recv directions.
+ # @debug
+ # def apply(self, variables):
+ # """Apply either mpi communications
+ # or local boundary conditions to fill ghosts.
+ # Loop over directions and switch among local BC or mpi comm.
+ # """
+ # assert isinstance(variables, list)
+ # local_bc_dir = self.topology.bc_dirs
+ # assert len(self.exchange_dir) + len(local_bc_dir) == self._dim
+
+ # i = 0
+ # for d in xrange(self._dim):
+ # if d in local_bc_dir:
+ # self._apply_bc_in_dir(variables, d)
+ # elif d in self.exchange_dir:
+ # self._apply_in_dir(variables, d, i)
+ # # update position in buffers list
+ # i += 1
+ # # End of loop through send/recv directions.
diff --git a/hysop/operator/discrete/particle_advection.py b/hysop/operator/discrete/particle_advection.py
index a8fb2be0945138cc3073879311e48c9aa5631266..d8cf10f58a5a6b4d3e0fc0d2bc0b13b68de07754 100644
--- a/hysop/operator/discrete/particle_advection.py
+++ b/hysop/operator/discrete/particle_advection.py
@@ -116,7 +116,7 @@ class ParticleAdvection(DiscreteOperator):
self._rhs_size, rwork=self._rw_integ,
f=num_interpolate,
topo=topo_fields,
- #optim=WITH_GUESS,
+ optim=WITH_GUESS,
indices=ic)
# -- remesh --
@@ -257,7 +257,7 @@ class ParticleAdvection(DiscreteOperator):
# - ghost-synchro of the rhs (if needed) is done by the odesolver
# - interpolation of velocity field from grid to particles is done
# 'inside' the ode solver (rhs = interpolation scheme)
- #self.num_advec.rwork[0][...] = self.velocity.data[self.direction][...]
+ self.num_advec.rwork[0][...] = self.velocity.data[self.direction][...]
self.part_positions = self.num_advec(
t, self.part_positions, dt, result=self.part_positions)