diff --git a/hysop/fields/ghost_exchangers.py b/hysop/fields/ghost_exchangers.py
index 6cc1d8a935b0dafbd47c03965d10890a2d18b614..e71f4e506dfd75234e563cf863d35bd480154464 100644
--- a/hysop/fields/ghost_exchangers.py
+++ b/hysop/fields/ghost_exchangers.py
@@ -1,4 +1,3 @@
-
from hysop.deps import np, hashlib
from hysop.tools.types import check_instance, to_tuple, first_not_None
from hysop.tools.numerics import default_invalid_value
@@ -18,7 +17,7 @@ class GhostExchanger(object):
check_instance(topology, TopologyView)
check_instance(data, tuple, minsize=1)
check_instance(ghost_op, GhostOperation)
-
+
dtype = data[0].dtype
for d in data[1:]:
if (d.dtype != dtype):
@@ -41,8 +40,8 @@ class GhostExchanger(object):
class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
- def __init__(self, name, topology, data,
- kind=None, directions=None,
+ def __init__(self, name, topology, data,
+ kind=None, directions=None,
ghosts=None, ghost_op=None,
exchange_method=None):
"""
@@ -53,26 +52,26 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
X P X
^
P < P > P
- v
+ v
X P X
-
+
Diagonal ghosts are set to NAN to ensure they are not used.
"""
ghost_op = first_not_None(ghost_op, GhostOperation.EXCHANGE)
-
- exchange_method = first_not_None(exchange_method,
+
+ exchange_method = first_not_None(exchange_method,
ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V)
- super(CartesianDiscreteFieldGhostExchanger, self).__init__(topology=topology,
+ super(CartesianDiscreteFieldGhostExchanger, self).__init__(topology=topology,
data=data, name=name, ghost_op=ghost_op)
-
+
mesh = self.mesh
dim = mesh.dim
directions = to_tuple(first_not_None(directions, range(dim)))
ghosts = to_tuple(first_not_None(ghosts, mesh.ghosts))
-
+
check_instance(topology, CartesianTopologyView)
check_instance(directions, tuple, values=int)
check_instance(ghosts, tuple, values=int, minval=0)
@@ -80,7 +79,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
assert len(set(directions)) == len(directions)
assert len(set(id(d) for d in data)) == len(data)
assert all(0<=d<dim for d in directions)
- assert len(ghosts)==dim or len(ghosts)==1
+ assert len(ghosts)==dim or len(ghosts)==1
if len(ghosts)==1:
ghosts*=dim
@@ -104,7 +103,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
msg='Unknown topology array backend kind {}.'.format(kind)
raise ValueError(msg)
self.kind = kind
-
+
self._launcher = None
def exchange_ghosts(self, **kwds):
@@ -135,14 +134,14 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
class _LauncherParameters(object):
def __init__(self):
# Call kwds depends on exchange method used (Send-Recv, Neighbor_AlltoAll, ...)
-
+
# SEND-RECV
# > one call per direction, per data buffer + (left + right)
self.isend_kwds = []
self.irecv_kwds = []
self.i_dst_buffers = []
self.i_src_buffers = []
-
+
# NEIGHBOR_ALL_TO_ALL_V
# > only one call for all buffers
self.v_kwds = None
@@ -152,7 +151,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
self.v_dst_buffers = ()
self.v_send_buffer = None
self.v_recv_buffer = None
-
+
# NEIGHBOR_ALL_TO_ALL_W
# > one call per data buffer
self.w_kwds = []
@@ -160,24 +159,24 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
self.w_dst_buffers = []
self.w_send_buffer = []
self.w_recv_buffer = []
-
+
# COMMON PARAMETERS
self.local_exchanges = []
self.diagonal_ghosts = []
self.has_mpi_exchanges = False
-
+
self.from_buffer = None # should source data be bufferized ?
self.to_buffer = None # should target data be bufferized ?
-
+
def msg_tag(i, local_rank, target_rank, d, direction):
tag = self.base_tag
tag += (i+1)*7919
- tag += (local_rank+1) * 967
+ tag += (local_rank+1) * 967
tag += (target_rank + 1) * 103
tag += (d + 1) * 97
tag += (direction+1)*17
return tag
-
+
dim = self.dim
ghosts = self.ghosts
ghost_op = self.ghost_op
@@ -191,7 +190,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
src_data_on_device = (self.kind is not Backend.HOST)
dst_data_on_device = (self.kind is not Backend.HOST)
-
+
all_outer_ghost_slices = self.all_outer_ghost_slices
lp = _LauncherParameters()
@@ -199,7 +198,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
lp.w_send_requests = {}
lp.w_recv_requests = {}
- if (ghost_op is GhostOperation.EXCHANGE) and (dim>1):
+ if (ghost_op is GhostOperation.EXCHANGE) and (dim>1) and not self.mesh.isExchInclDiag:
# set all outer diagonal ghosts to NAN
for directions in all_outer_ghost_slices[dim]:
if all(ghosts[d]>0 for d in directions):
@@ -216,10 +215,10 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
continue
inner_left, inner_right, shape = self.inner_ghosts[d]
outer_left, outer_right, shape = self.outer_ghosts[d]
-
+
left_rank = neighbour_ranks[0,d]
right_rank = neighbour_ranks[1,d]
-
+
nprocs = proc_shape[d]
lp.has_mpi_exchanges |= (nprocs > 1)
@@ -233,7 +232,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
# OPERATION IS
# OUTER_GHOSTS[...] = INNER_GHOSTS
if (exchange_method is ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V):
- # Send and receive every buffer at once to neighbours
+ # Send and receive every buffer at once to neighbours
# /!\ we send and receive all data components at once
if nprocs==2:
# switch left and right in 2 proc periodic case
@@ -267,7 +266,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
elif (exchange_method is ExchangeMethod.ISEND_IRECV):
# Exchanges with left neighour
assert (left_rank != local_rank) and (left_rank != -1), left_rank
-
+
# send inner left to left rank
sendtag = msg_tag(i, local_rank, left_rank, d, 0)
if src_data_on_device:
@@ -280,11 +279,11 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
send_buf = buf.handle
mpi_type = TopoTools.create_subarray(inner_left, buf.shape,
mpi_type=base_mpi_type)
- send_kwds = {'buf':[send_buf, mpi_type],
- 'dest':left_rank,
+ send_kwds = {'buf':[send_buf, mpi_type],
+ 'dest':left_rank,
'tag':sendtag}
lp.isend_kwds.append(send_kwds)
-
+
# receive outer right from left rank
recvtag = msg_tag(i, left_rank, local_rank, d, 1)
if dst_data_on_device:
@@ -297,11 +296,11 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
recv_buf = buf.handle
mpi_type = TopoTools.create_subarray(outer_left, buf.shape,
mpi_type=base_mpi_type)
- recv_kwds = {'buf':[recv_buf, mpi_type],
- 'source':left_rank,
+ recv_kwds = {'buf':[recv_buf, mpi_type],
+ 'source':left_rank,
'tag':recvtag}
lp.irecv_kwds.append(recv_kwds)
-
+
# Exchanges with right neighour
assert (right_rank != local_rank) and (right_rank != -1)
@@ -317,11 +316,11 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
send_buf = buf.handle
mpi_type = TopoTools.create_subarray(inner_right, buf.shape,
mpi_type=base_mpi_type)
- send_kwds = {'buf':[send_buf, mpi_type],
- 'dest':right_rank,
+ send_kwds = {'buf':[send_buf, mpi_type],
+ 'dest':right_rank,
'tag':sendtag}
lp.isend_kwds.append(send_kwds)
-
+
# receive outer left from right rank
recvtag = msg_tag(i, right_rank, local_rank, d, 0)
if dst_data_on_device:
@@ -334,8 +333,8 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
recv_buf = buf.handle
mpi_type = TopoTools.create_subarray(outer_right, buf.shape,
mpi_type=base_mpi_type)
- recv_kwds = {'buf':[recv_buf, mpi_type],
- 'source':right_rank,
+ recv_kwds = {'buf':[recv_buf, mpi_type],
+ 'source':right_rank,
'tag':recvtag}
lp.irecv_kwds.append(recv_kwds)
lp.from_buffer = src_data_on_device
@@ -349,7 +348,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
# OPERATION IS
# INNER GHOSTS = OP(INNER_GHOSTS, OUTER_GHOSTS)
if (exchange_method is ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V):
- # Send and receive every buffer at once to neighbours
+ # Send and receive every buffer at once to neighbours
# /!\ we send and receive all data components at once
left_rank = neighbour_ranks[0,d]
right_rank = neighbour_ranks[1,d]
@@ -388,7 +387,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
# Exchanges with left neighour
left_rank = neighbour_ranks[0,d]
assert (left_rank != local_rank) and (left_rank != -1)
-
+
# send outer left to left rank
sendtag = msg_tag(i, local_rank, left_rank, d, 0)
if src_data_on_device:
@@ -401,26 +400,26 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
send_buf = buf.handle
mpi_type = TopoTools.create_subarray(outer_left, buf.shape,
mpi_type=base_mpi_type)
- send_kwds = {'buf':[send_buf, mpi_type],
- 'dest':left_rank,
+ send_kwds = {'buf':[send_buf, mpi_type],
+ 'dest':left_rank,
'tag':sendtag}
lp.isend_kwds.append(send_kwds)
-
+
# receive outer right ghosts data from left rank in a tmp buffer
recvtag = msg_tag(i, left_rank, local_rank, d, 1)
- mpi_type = base_mpi_type.Create_contiguous(buf[inner_left].size)
+ mpi_type = base_mpi_type.Create_contiguous(buf[inner_left].size)
mpi_type.Commit()
tmp = self.host_backend.empty(shape=shape, dtype=base_dtype)
- recv_kwds = {'buf':[tmp.handle, mpi_type],
- 'source':left_rank,
+ recv_kwds = {'buf':[tmp.handle, mpi_type],
+ 'source':left_rank,
'tag':recvtag}
lp.irecv_kwds.append(recv_kwds)
lp.i_dst_buffers.append((tmp, buf, inner_left))
-
+
# Exchanges with right neighour
right_rank = neighbour_ranks[1,d]
assert (right_rank != local_rank) and (right_rank != -1)
-
+
# send outer right to right rank
sendtag = msg_tag(i, local_rank, right_rank, d, 1)
if src_data_on_device:
@@ -433,18 +432,18 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
send_buf = buf.handle
mpi_type = TopoTools.create_subarray(outer_right, buf.shape,
mpi_type=base_mpi_type)
- send_kwds = {'buf':[send_buf, mpi_type],
- 'dest':right_rank,
+ send_kwds = {'buf':[send_buf, mpi_type],
+ 'dest':right_rank,
'tag':sendtag}
lp.isend_kwds.append(send_kwds)
-
+
# receive outer left ghosts data from right rank in a tmp buffer
recvtag = msg_tag(i, right_rank, local_rank, d, 0)
- mpi_type = base_mpi_type.Create_contiguous(buf[inner_right].size)
+ mpi_type = base_mpi_type.Create_contiguous(buf[inner_right].size)
mpi_type.Commit()
tmp = self.host_backend.empty(shape=shape, dtype=base_dtype)
- recv_kwds = {'buf':[tmp.handle, mpi_type],
- 'source':right_rank,
+ recv_kwds = {'buf':[tmp.handle, mpi_type],
+ 'source':right_rank,
'tag':recvtag}
lp.irecv_kwds.append(recv_kwds)
lp.i_dst_buffers.append((tmp, buf, inner_right))
@@ -464,7 +463,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
sdispls, rdispls = (),() # in bytes
sdispl, rdispl = 0, 0
scount, rcount = 0, 0
-
+
src_buffers, dst_buffers = (), ()
for direction in xrange(comm.Get_dim()):
for (tag,rank) in zip(('left','right'), comm.Shift(direction, 1)):
@@ -482,7 +481,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vslc = slice(scount, scount+send_count)
src_buffers += ((src[slc],vslc),)
else:
- mpi_type = TopoTools.create_subarray(slc, src.shape,
+ mpi_type = TopoTools.create_subarray(slc, src.shape,
mpi_type=base_mpi_type)
send_count = 1
else:
@@ -493,7 +492,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
sdispls += (sdispl,)
sdispl += displ
scount += send_count
-
+
displ = 0
if ((tag,rank) in lp.w_recv_requests):
requests = lp.w_recv_requests[(tag,rank)]
@@ -508,7 +507,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vslc = slice(rcount, rcount+recv_count)
dst_buffers += ((dst[slc],vslc),)
else:
- mpi_type = TopoTools.create_subarray(slc, dst.shape,
+ mpi_type = TopoTools.create_subarray(slc, dst.shape,
mpi_type=base_mpi_type)
recv_count = 1
else:
@@ -543,12 +542,12 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
'sendbuf': [sendbuf, sendcounts, sdispls, sendtypes],
'recvbuf': [recvbuf, recvcounts, rdispls, recvtypes]
} )
-
+
# handle all_to_all_v kwds (one call for all buffers and all neighbours)
if lp.has_mpi_exchanges and (exchange_method is ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V):
assert lp.v_send_requests
assert lp.v_recv_requests
- sendtype, recvtype = base_mpi_type, base_mpi_type
+ sendtype, recvtype = base_mpi_type, base_mpi_type
sendcounts, recvcounts = (),() # in type counts
sdispls, rdispls = (),() # in type counts
send_disp, recv_disp = 0,0 # in type counts
@@ -591,18 +590,18 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
lp.v_recv_buffer = recv_buffer
return lp
-
+
def _build_python_launcher(self):
assert (self._launcher is None)
-
+
ghost_op = self.ghost_op
exchange_method = self.exchange_method
comm = self.topology.comm
base_dtype = self.base_dtype
lp = self._prepare_launcher()
-
+
if ghost_op is GhostOperation.EXCHANGE:
def local_exchanges():
for (buf,outer,inner,shape,direction) in lp.local_exchanges:
@@ -716,7 +715,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
raise NotImplementedError(msg)
self._launcher = python_launcher
-
+
def _build_opencl_launcher(self):
assert (self._launcher is None)
@@ -729,14 +728,14 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
self._fn = fn
def __call__(self, *args, **kwds):
return self._fn(*args, **kwds)
-
+
lp = self._prepare_launcher()
dim = self.dim
ghost_op = self.ghost_op
exchange_method = self.exchange_method
comm = self.topology.comm
base_dtype = self.base_dtype
-
+
# generate the minimal number of temporary backend buffers
tmp_buffers = {}
def mk_tmp(shape, dtype, color=0):
@@ -750,12 +749,12 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
return tmp
if ghost_op is GhostOperation.EXCHANGE:
- local_kl = OpenClKernelListLauncher(name='local_ghosts_exchanges_{}'.format(self.name))
+ local_kl = OpenClKernelListLauncher(name='local_ghosts_exchanges_{}'.format(self.name))
for i,(buf,outer_slc,inner_slc,shape,direction) in enumerate(lp.local_exchanges):
dirlabel = DirectionLabels[dim-direction-1]
vname = '{}_{}_{}'.format(self.name, i, dirlabel)
-
- # some opencl platforms reject inplace buffer copies
+
+ # some opencl platforms reject inplace buffer copies
# so we use a tmp buffer
# to perform local ghost exchanges
tmp = mk_tmp(shape=shape, dtype=base_dtype)
@@ -765,7 +764,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
varname=vname+'_inner_to_tmp_ghost_buffer',
src=buf, src_slices=inner_slc, dst=tmp)
local_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_tmp_ghost_buffer_to_outer',
+ varname=vname+'_tmp_ghost_buffer_to_outer',
src=tmp, dst=buf, dst_slices=outer_slc)
for i,(buf,slc,shape,value) in enumerate(lp.diagonal_ghosts):
@@ -774,8 +773,8 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
tmp.fill(value)
local_kl += OpenClCopyBufferRectLauncher.from_slices(
varname=vname, src=tmp, dst=buf, dst_slices=slc)
-
- kl = OpenClKernelListLauncher(name='exchange_ghosts_{}'.format(self.name))
+
+ kl = OpenClKernelListLauncher(name='exchange_ghosts_{}'.format(self.name))
if lp.has_mpi_exchanges:
assert (lp.from_buffer is True)
assert (lp.to_buffer is True)
@@ -783,7 +782,7 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
assert len(lp.w_kwds)>0
d2h_kls = []
h2d_kls = []
-
+
for i,(send_buffer,src_buffers) in \
enumerate(zip(lp.w_send_buffer, lp.w_src_buffers)):
d2h_kl_i = OpenClKernelListLauncher(
@@ -792,11 +791,11 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vname = '{}_{}_{}'.format(self.name, i, j)
dst = send_buffer[dst_slc].reshape(src.shape)
d2h_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_ghosts_D2H',
- src=src, src_slices=None,
+ varname=vname+'_inner_ghosts_D2H',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
d2h_kls.append(d2h_kl_i)
-
+
for i,(recv_buffer,dst_buffers) in \
enumerate(zip(lp.w_recv_buffer, lp.w_dst_buffers)):
h2d_kl_i = OpenClKernelListLauncher(
@@ -804,8 +803,8 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
for j,(dst,src_slc) in enumerate(dst_buffers):
src = recv_buffer[src_slc].reshape(dst.shape)
h2d_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_ghosts_H2D',
- src=src, src_slices=None,
+ varname=vname+'_outer_ghosts_H2D',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
h2d_kls.append(h2d_kl_i)
@@ -833,18 +832,18 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vname = '{}_{}'.format(self.name, i)
dst = lp.v_send_buffer[dst_slc].reshape(src.shape)
d2h_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_ghosts_D2H',
- src=src, src_slices=None,
+ varname=vname+'_inner_ghosts_D2H',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
-
+
for i,(dst,src_slc) in enumerate(lp.v_dst_buffers):
vname = '{}_{}'.format(self.name, i)
src = lp.v_recv_buffer[src_slc].reshape(dst.shape)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_ghosts_H2D',
- src=src, src_slices=None,
+ varname=vname+'_outer_ghosts_H2D',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
-
+
def opencl_launcher(**kwds):
evt = d2h_kl(**kwds)
evt.wait()
@@ -860,20 +859,20 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
assert len(lp.isend_kwds) == len(lp.i_src_buffers)
d2h_kls = []
h2d_kls = []
-
+
for i,(dst, src, src_slc) in enumerate(lp.i_src_buffers):
vname = '{}_{}'.format(self.name, i)
d2h_kl = OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_ghosts_D2H',
- src=src, src_slices=src_slc,
+ varname=vname+'_inner_ghosts_D2H',
+ src=src, src_slices=src_slc,
dst=dst, dst_slices=None)
d2h_kls.append(d2h_kl)
-
+
for i,(src, dst, dst_slc) in enumerate(lp.i_dst_buffers):
vname = '{}_{}'.format(self.name, i)
h2d_kl = OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_ghosts_H2D',
- src=src, src_slices=None,
+ varname=vname+'_outer_ghosts_H2D',
+ src=src, src_slices=None,
dst=dst, dst_slices=dst_slc)
h2d_kls.append(h2d_kl)
@@ -901,34 +900,34 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
else:
kl += local_kl
elif ghost_op is GhostOperation.ACCUMULATE:
- local_kl = OpenClKernelListLauncher(name='local_ghosts_accumulate_{}'.format(self.name))
+ local_kl = OpenClKernelListLauncher(name='local_ghosts_accumulate_{}'.format(self.name))
for i,(buf,outer_slc,inner_slc,shape,direction) in enumerate(lp.local_exchanges):
dirlabel = DirectionLabels[dim-direction-1]
vname = '{}_{}_{}'.format(self.name, i, dirlabel)
-
+
ltmp = mk_tmp(shape=shape, dtype=base_dtype, color=0)
rtmp = mk_tmp(shape=shape, dtype=base_dtype, color=1)
# accumulate outer ghosts to inner ghosts
local_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_to_tmp_ghost_buffer_0',
+ varname=vname+'_inner_to_tmp_ghost_buffer_0',
src=buf, src_slices=inner_slc, dst=ltmp)
local_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_to_tmp_ghost_buffer_1',
+ varname=vname+'_outer_to_tmp_ghost_buffer_1',
src=buf, src_slices=outer_slc, dst=rtmp)
- local_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
+ local_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
build_kernel_launcher=True)
local_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_accumulation_to_inner',
+ varname=vname+'_accumulation_to_inner',
src=ltmp, dst=buf, dst_slices=inner_slc)
-
- kl = OpenClKernelListLauncher(name='accumulate_ghosts_{}'.format(self.name))
+
+ kl = OpenClKernelListLauncher(name='accumulate_ghosts_{}'.format(self.name))
if lp.has_mpi_exchanges:
if exchange_method is ExchangeMethod.NEIGHBOR_ALL_TO_ALL_W:
assert len(lp.w_kwds)>0
d2h_kls = []
h2d_kls = []
-
+
for i,(send_buffer,src_buffers) in \
enumerate(zip(lp.w_send_buffer, lp.w_src_buffers)):
d2h_kl_i = OpenClKernelListLauncher(
@@ -937,11 +936,11 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vname = '{}_{}_{}'.format(self.name, i, j)
dst = send_buffer[dst_slc].reshape(src.shape)
d2h_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_ghosts_D2H',
- src=src, src_slices=None,
+ varname=vname+'_inner_ghosts_D2H',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
d2h_kls.append(d2h_kl_i)
-
+
for i,(recv_buffer,dst_buffers) in \
enumerate(zip(lp.w_recv_buffer, lp.w_dst_buffers)):
h2d_kl_i = OpenClKernelListLauncher(
@@ -954,17 +953,17 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
ltmp = mk_tmp(shape=shape, dtype=base_dtype, color=0)
rtmp = mk_tmp(shape=shape, dtype=base_dtype, color=1)
h2d_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_to_tmp_ghost_buffer_0',
- src=outer, src_slices=None,
+ varname=vname+'_outer_to_tmp_ghost_buffer_0',
+ src=outer, src_slices=None,
dst=ltmp, dst_slices=None)
h2d_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_to_tmp_ghost_buffer_1',
- src=inner, src_slices=None,
+ varname=vname+'_inner_to_tmp_ghost_buffer_1',
+ src=inner, src_slices=None,
dst=rtmp, dst_slices=None)
- h2d_kl_i += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
+ h2d_kl_i += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
build_kernel_launcher=True)
h2d_kl_i += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_accumulation_to_inner',
+ varname=vname+'_accumulation_to_inner',
src=ltmp, dst=inner)
h2d_kls.append(h2d_kl_i)
@@ -994,10 +993,10 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
vname = '{}_{}'.format(self.name, i)
dst = lp.v_send_buffer[dst_slc].reshape(src.shape)
d2h_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_ghosts_D2H',
- src=src, src_slices=None,
+ varname=vname+'_inner_ghosts_D2H',
+ src=src, src_slices=None,
dst=dst, dst_slices=None)
-
+
for i,(dst,src_slc) in enumerate(lp.v_dst_buffers):
vname = '{}_{}'.format(self.name, i)
shape = dst.shape
@@ -1006,17 +1005,17 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
ltmp = mk_tmp(shape=shape, dtype=base_dtype, color=0)
rtmp = mk_tmp(shape=shape, dtype=base_dtype, color=1)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_to_tmp_ghost_buffer_0',
- src=outer, src_slices=None,
+ varname=vname+'_outer_to_tmp_ghost_buffer_0',
+ src=outer, src_slices=None,
dst=ltmp, dst_slices=None)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_to_tmp_ghost_buffer_1',
- src=inner, src_slices=None,
+ varname=vname+'_inner_to_tmp_ghost_buffer_1',
+ src=inner, src_slices=None,
dst=rtmp, dst_slices=None)
- h2d_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
+ h2d_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
build_kernel_launcher=True)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_accumulation_to_inner',
+ varname=vname+'_accumulation_to_inner',
src=ltmp, dst=inner)
def opencl_launcher(**kwds):
@@ -1034,15 +1033,15 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
assert len(lp.isend_kwds) == len(lp.i_src_buffers)
d2h_kls = []
h2d_kls = []
-
+
for i,(dst, src, src_slc) in enumerate(lp.i_src_buffers):
vname = '{}_{}'.format(self.name, i)
d2h_kl = OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_ghosts_D2H',
- src=src, src_slices=src_slc,
+ varname=vname+'_outer_ghosts_D2H',
+ src=src, src_slices=src_slc,
dst=dst, dst_slices=None)
d2h_kls.append(d2h_kl)
-
+
for i,(src, dst, dst_slc) in enumerate(lp.i_dst_buffers):
shape = src.shape
outer = src
@@ -1053,17 +1052,17 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
h2d_kl = OpenClKernelListLauncher(
'outer_ghosts_H2D_exchanges_{}_{}'.format(self.name,i))
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_outer_to_tmp_ghost_buffer_0',
- src=outer, src_slices=None,
+ varname=vname+'_outer_to_tmp_ghost_buffer_0',
+ src=outer, src_slices=None,
dst=ltmp, dst_slices=None)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_inner_to_tmp_ghost_buffer_1',
- src=inner, src_slices=None,
+ varname=vname+'_inner_to_tmp_ghost_buffer_1',
+ src=inner, src_slices=None,
dst=rtmp, dst_slices=None)
- h2d_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
+ h2d_kl += self.backend.add(x1=ltmp, x2=rtmp, out=ltmp,
build_kernel_launcher=True)
h2d_kl += OpenClCopyBufferRectLauncher.from_slices(
- varname=vname+'_accumulation_to_inner',
+ varname=vname+'_accumulation_to_inner',
src=ltmp, dst=inner)
h2d_kls.append(h2d_kl)
@@ -1089,7 +1088,6 @@ class CartesianDiscreteFieldGhostExchanger(GhostExchanger):
msg='Unknown MPI exchange method {}.'.format(exchange_method)
raise NotImplementedError(msg)
else:
- kl += local_kl
+ kl += local_kl
self._launcher = kl
-
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index 507b4bcea01d918b5f515e70d464f48ddd197f36..16e2d8ae85d13389eee50d6e705537086b86db41 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -10,16 +10,16 @@ from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_topology import CartesianTopology, CartesianTopologyState
from hysop.testsenv import iter_clenv, test_context
from hysop.tools.numerics import is_fp, is_integer
-
+
def __random_init(data, coords):
shape = data[0].shape
dtype = data[0].dtype
if is_integer(dtype):
for d in data:
- d[...] = npw.random.random_integers(low=0, high=255, size=shape)
+ d[...] = npw.random.random_integers(low=0, high=255, size=shape)
elif is_fp(dtype):
for d in data:
- d[...] = npw.random.random(size=d.shape)
+ d[...] = npw.random.random(size=d.shape)
else:
msg = 'Unknown dtype {}.'.format(dtype)
raise NotImplementedError(msg)
@@ -40,16 +40,16 @@ def test_serial_initialization_1d():
npts = (10,)
nghosts = (2,)
discretization = Discretization(npts, nghosts)
-
+
domain = Box(dim=1)
- F0 = Field(domain=domain, name='F0', nb_components=1)
- F1 = Field(domain=domain, name='F1', nb_components=2)
-
- topo0 = CartesianTopology(domain=domain, discretization=discretization,
+ F0 = Field(domain=domain, name='F0', nb_components=1)
+ F1 = Field(domain=domain, name='F1', nb_components=2)
+
+ topo0 = CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.HOST)
topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
-
+
for topo in topos:
print ' Topo {}::{} | {}'.format(topo.full_tag, topo.backend.kind, topo.backend)
dF0 = F0.discretize(topo)
@@ -71,16 +71,16 @@ def test_serial_initialization_2d():
npts = (11,13)
nghosts = (3,5)
discretization = Discretization(npts, nghosts)
-
+
domain = Box(dim=2)
- F0 = Field(domain=domain, name='F0', nb_components=1)
- F1 = Field(domain=domain, name='F1', nb_components=2)
-
- topo0 = CartesianTopology(domain=domain, discretization=discretization,
+ F0 = Field(domain=domain, name='F0', nb_components=1)
+ F1 = Field(domain=domain, name='F1', nb_components=2)
+
+ topo0 = CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.HOST)
topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
-
+
for topo in topos:
print ' Topo {}::{} | {}'.format(topo.full_tag, topo.backend.kind, topo.backend)
dF0 = F0.discretize(topo)
@@ -92,11 +92,39 @@ def test_serial_initialization_2d():
for (i,d) in enumerate(data):
print ' *buffer {}'.format(i)
buf = d.get().handle
+ slo_x = (slice(0,nghosts[0]),
+ slice(nghosts[0],nghosts[0]+npts[0]-1),
+ slice(nghosts[0]+npts[0]-1,None))
+ slo_y = (slice(0,nghosts[1]),
+ slice(nghosts[1],nghosts[1]+npts[1]-1),
+ slice(nghosts[1]+npts[1]-1,None))
+ sli_x = (slice(nghosts[0], 2*nghosts[0]),
+ slice(2*nghosts[0],npts[0]-1),
+ slice(npts[0]-1, nghosts[0]+npts[0]-1))
+ sli_y = (slice(nghosts[1], 2*nghosts[1]),
+ slice(2*nghosts[1],npts[1]-1),
+ slice(npts[1]-1, nghosts[1]+npts[1]-1))
assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1)
- assert (buf[nghosts[0]:2*nghosts[0],:] == buf[nghosts[0]+npts[0]-1:,:]).all()
- assert (buf[:nghosts[0],:] == buf[npts[0]-1:nghosts[0]+npts[0]-1,:]).all()
- assert (buf[:,nghosts[1]:2*nghosts[1]] == buf[:,nghosts[1]+npts[1]-1:]).all()
- assert (buf[:,:nghosts[1]] == buf[:,npts[1]-1:nghosts[1]+npts[1]-1]).all()
+ assert buf[slo_x[0],slo_y[0]].shape == nghosts
+ assert buf[slo_x[2],slo_y[2]].shape == nghosts
+ assert buf[slo_x[1],slo_y[1]].shape == tuple([_-1 for _ in npts])
+ assert buf[sli_x[0],sli_y[0]].shape == nghosts
+ assert buf[sli_x[2],sli_y[2]].shape == nghosts
+ assert np.all(buf[slo_x[0],slo_y[1]] == buf[sli_x[2],slo_y[1]])
+ assert np.all(buf[slo_x[2],slo_y[1]] == buf[sli_x[0],slo_y[1]])
+ assert np.all(buf[slo_x[1],slo_y[0]] == buf[slo_x[1],sli_y[2]])
+ assert np.all(buf[slo_x[1],slo_y[2]] == buf[slo_x[1],sli_y[0]])
+ if topo.mesh.isExchInclDiag:
+ assert np.all(buf[slo_x[0],slo_y[0]]==buf[sli_x[2],sli_y[2]])
+ assert np.all(buf[slo_x[2],slo_y[0]]==buf[sli_x[0],sli_y[2]])
+ assert np.all(buf[slo_x[2],slo_y[2]]==buf[sli_x[0],sli_y[0]])
+ assert np.all(buf[slo_x[0],slo_y[2]]==buf[sli_x[2],sli_y[0]])
+ else:
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[0]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[0]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[2]]))
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[2]]))
+
def test_serial_initialization_3d():
print
@@ -104,16 +132,16 @@ def test_serial_initialization_3d():
npts = (17,11,15)
nghosts = (5,2,3)
discretization = Discretization(npts, nghosts)
-
+
domain = Box(dim=3)
- F0 = Field(domain=domain, name='F0', nb_components=1)
- F1 = Field(domain=domain, name='F1', nb_components=2)
-
- topo0 = CartesianTopology(domain=domain, discretization=discretization,
+ F0 = Field(domain=domain, name='F0', nb_components=1)
+ F1 = Field(domain=domain, name='F1', nb_components=2)
+
+ topo0 = CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.HOST)
topos = (topo0,) + tuple(CartesianTopology(domain=domain, discretization=discretization,
backend=Backend.OPENCL, cl_env=cl_env) for cl_env in iter_clenv())
-
+
for topo in topos:
print ' Topo {}::{} | {}'.format(topo.full_tag, topo.backend.kind, topo.backend)
dF0 = F0.discretize(topo)
@@ -125,14 +153,60 @@ def test_serial_initialization_3d():
for (i,d) in enumerate(data):
print ' *buffer {}'.format(i)
buf = d.get().handle
+ slo_x = (slice(0,nghosts[0]),
+ slice(nghosts[0],nghosts[0]+npts[0]-1),
+ slice(nghosts[0]+npts[0]-1,None))
+ slo_y = (slice(0,nghosts[1]),
+ slice(nghosts[1],nghosts[1]+npts[1]-1),
+ slice(nghosts[1]+npts[1]-1,None))
+ slo_z = (slice(0,nghosts[2]),
+ slice(nghosts[2],nghosts[2]+npts[2]-1),
+ slice(nghosts[2]+npts[2]-1,None))
+ sli_x = (slice(nghosts[0], 2*nghosts[0]),
+ slice(2*nghosts[0],npts[0]-1),
+ slice(npts[0]-1, nghosts[0]+npts[0]-1))
+ sli_y = (slice(nghosts[1], 2*nghosts[1]),
+ slice(2*nghosts[1],npts[1]-1),
+ slice(npts[1]-1, nghosts[1]+npts[1]-1))
+ sli_z = (slice(nghosts[2], 2*nghosts[2]),
+ slice(2*nghosts[2],npts[2]-1),
+ slice(npts[2]-1, nghosts[2]+npts[2]-1))
assert buf.shape==(npts[0]+2*nghosts[0]-1,npts[1]+2*nghosts[1]-1,npts[2]
+2*nghosts[2]-1)
- assert (buf[nghosts[0]:2*nghosts[0],:,:] == buf[nghosts[0]+npts[0]-1:,:,:]).all()
- assert (buf[:nghosts[0],:,:] == buf[npts[0]-1:nghosts[0]+npts[0]-1,:,:]).all()
- assert (buf[:,nghosts[1]:2*nghosts[1],:] == buf[:,nghosts[1]+npts[1]-1:,:]).all()
- assert (buf[:,:nghosts[1],:] == buf[:,npts[1]-1:nghosts[1]+npts[1]-1,:]).all()
- assert (buf[:,:,nghosts[2]:2*nghosts[2]] == buf[:,:,nghosts[2]+npts[2]-1:]).all()
- assert (buf[:,:,:nghosts[2]] == buf[:,:,npts[2]-1:nghosts[2]+npts[2]-1]).all()
+ assert buf[slo_x[0],slo_y[0],slo_z[0]].shape == nghosts
+ assert buf[slo_x[2],slo_y[2],slo_z[0]].shape == nghosts
+ assert buf[slo_x[0],slo_y[0],slo_z[2]].shape == nghosts
+ assert buf[slo_x[2],slo_y[2],slo_z[2]].shape == nghosts
+ assert buf[slo_x[1],slo_y[1],slo_z[1]].shape == tuple([_-1 for _ in npts])
+ assert buf[sli_x[0],sli_y[0],sli_z[0]].shape == nghosts
+ assert buf[sli_x[2],sli_y[2],sli_z[0]].shape == nghosts
+ assert buf[sli_x[0],sli_y[0],sli_z[2]].shape == nghosts
+ assert buf[sli_x[2],sli_y[2],sli_z[2]].shape == nghosts
+ assert np.all(buf[slo_x[0],slo_y[1],slo_z[1]] == buf[sli_x[2],slo_y[1],slo_z[1]])
+ assert np.all(buf[slo_x[2],slo_y[1],slo_z[1]] == buf[sli_x[0],slo_y[1],slo_z[1]])
+ assert np.all(buf[slo_x[1],slo_y[0],slo_z[1]] == buf[slo_x[1],sli_y[2],slo_z[1]])
+ assert np.all(buf[slo_x[1],slo_y[2],slo_z[1]] == buf[slo_x[1],sli_y[0],slo_z[1]])
+ assert np.all(buf[slo_x[1],slo_y[1],slo_z[0]] == buf[slo_x[1],slo_y[1],sli_z[2]])
+ assert np.all(buf[slo_x[1],slo_y[1],slo_z[2]] == buf[slo_x[1],slo_y[1],sli_z[0]])
+ if topo.mesh.isExchInclDiag:
+ assert np.all(buf[slo_x[0],slo_y[0],slo_z[0]]==buf[sli_x[2],sli_y[2],sli_z[2]])
+ assert np.all(buf[slo_x[2],slo_y[0],slo_z[0]]==buf[sli_x[0],sli_y[2],sli_z[2]])
+ assert np.all(buf[slo_x[2],slo_y[2],slo_z[0]]==buf[sli_x[0],sli_y[0],sli_z[2]])
+ assert np.all(buf[slo_x[0],slo_y[2],slo_z[0]]==buf[sli_x[2],sli_y[0],sli_z[2]])
+ assert np.all(buf[slo_x[0],slo_y[0],slo_z[2]]==buf[sli_x[2],sli_y[2],sli_z[0]])
+ assert np.all(buf[slo_x[2],slo_y[0],slo_z[2]]==buf[sli_x[0],sli_y[2],sli_z[0]])
+ assert np.all(buf[slo_x[2],slo_y[2],slo_z[2]]==buf[sli_x[0],sli_y[0],sli_z[0]])
+ assert np.all(buf[slo_x[0],slo_y[2],slo_z[2]]==buf[sli_x[2],sli_y[0],sli_z[0]])
+ else:
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[0],slo_z[0]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[0],slo_z[0]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[2],slo_z[0]]))
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[2],slo_z[0]]))
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[0],slo_z[2]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[0],slo_z[2]]))
+ assert np.all(np.isnan(buf[slo_x[2],slo_y[2],slo_z[2]]))
+ assert np.all(np.isnan(buf[slo_x[0],slo_y[2],slo_z[2]]))
+
def iter_backends():
yield (Backend.HOST, None)
@@ -157,35 +231,35 @@ def test_mpi_ghost_exchange(comm):
for dim in xrange(1,4+__ENABLE_LONG_TESTS__):
if rank==0:
print(' >DIM={}'.format(dim))
-
+
npts = (53,47,59,23)[:dim]
nghosts = (2,1,0,3)[:dim]
discretization = Discretization(npts, nghosts)
domain = Box(dim=dim)
-
+
for dtype in dtypes[size-1:size]:
if rank==0:
print(' >DTYPE={}'.format(dtype))
-
- F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
- F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
-
+
+ F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
+ F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
+
for (backend, cl_env) in iter_backends():
if rank==0:
- print ' >BACKEND.{}{}'.format(backend,
+ print ' >BACKEND.{}{}'.format(backend,
'' if (cl_env is None) else '::{}.{}'.format(
cl_env.platform.name.strip(), cl_env.device.name.strip()))
- for shape in it.product(xrange(0,size+1), repeat=dim):
+ for shape in it.product(xrange(0,size+1), repeat=dim):
if np.prod(shape, dtype=np.uint32)!=size:
continue
if rank==0:
print ' *cart_shape: {}'.format(shape)
- topo = CartesianTopology(domain=domain, discretization=discretization,
+ topo = CartesianTopology(domain=domain, discretization=discretization,
backend=backend, cart_shape=shape, cl_env=cl_env)
assert (topo.proc_shape==shape).all()
def ghost_base(i,d,rank,local_dir):
- return (i+1)*17 + (d+1)*13 + (local_dir+1)*11 + (rank+1)*7
+ return (i+1)*17 + (d+1)*13 + (local_dir+1)*11 + (rank+1)*7
def ghost_vals(shape,dtype,i,d,rank,local_dir):
if (shape is None) or (len(shape)==0):
raise ValueError('Shape is None or an empty tuple: {}'.format(shape))
@@ -193,7 +267,7 @@ def test_mpi_ghost_exchange(comm):
I = np.ix_(*tuple(np.arange(shape[direction], dtype=dtype) \
for direction in xrange(len(shape))))
return base + I[d]
-
+
for F in (F0, F1):
dF = F.discretize(topo)
if rank==0:
@@ -207,19 +281,19 @@ def test_mpi_ghost_exchange(comm):
for d in xrange(dim):
dF.initialize(__zero_init, exchange_ghosts=False)
if rank==0:
- print DirectionLabels[dim-1-d],
+ print DirectionLabels[dim-1-d],
sys.stdout.flush()
lghosts, rghosts, shape = dF.inner_ghost_slices[d]
_lghosts, _rghosts, shape = dF.outer_ghost_slices[d]
for (i,data) in enumerate(dF.data):
data[lghosts] = ghost_vals(shape, dtype, i,d,rank,0)
- data[rghosts] = ghost_vals(shape, dtype, i,d,rank,1)
+ data[rghosts] = ghost_vals(shape, dtype, i,d,rank,1)
data[_lghosts] = -10
data[_rghosts] = +10
-
+
dF.exchange_ghosts(directions=d, exchange_method=exchange_method)
-
+
lghosts, rghosts, shape = dF.outer_ghost_slices[d]
left_rank, right_rank = topo.proc_neighbour_ranks[:,d]
if (left_rank == -1):
@@ -251,7 +325,7 @@ def test_mpi_ghost_accumulate(comm):
print msg
print '*'*len(msg)
print 'test_mpi_ghost_accumulate()'.format(size)
-
+
dtypes = (np.float32, np.float32, np.float64,
np.complex64, np.complex128,
np.int16, np.int32, np.int64,
@@ -260,30 +334,30 @@ def test_mpi_ghost_accumulate(comm):
for dim in xrange(1,4+__ENABLE_LONG_TESTS__):
if rank==0:
print(' >DIM={}'.format(dim))
-
+
npts = (53,57,51,49)[:dim]
nghosts = (1,3,0,2)[:dim]
discretization = Discretization(npts, nghosts)
domain = Box(dim=dim)
-
+
for dtype in dtypes[size-1:size]:
if rank==0:
print(' >DTYPE={}'.format(dtype))
-
- F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
- F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
-
+
+ F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
+ F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
+
for (backend, cl_env) in iter_backends():
if rank==0:
- print ' >BACKEND.{}{}'.format(backend,
+ print ' >BACKEND.{}{}'.format(backend,
'' if (cl_env is None) else '::{}.{}'.format(
cl_env.platform.name.strip(), cl_env.device.name.strip()))
- for shape in it.product(xrange(0,size+1), repeat=dim):
+ for shape in it.product(xrange(0,size+1), repeat=dim):
if np.prod(shape, dtype=np.uint32)!=size:
continue
if rank==0:
print ' *cart_shape: {}'.format(shape)
- topo = CartesianTopology(domain=domain, discretization=discretization,
+ topo = CartesianTopology(domain=domain, discretization=discretization,
backend=backend, cart_shape=shape, cl_env=cl_env)
assert (topo.proc_shape==shape).all()
@@ -292,7 +366,7 @@ def test_mpi_ghost_accumulate(comm):
dirweight = np.asarray([37,41,43,51], dtype=np.int32)
directions = np.asarray(directions, dtype=np.int32) + 1
displacements = np.asarray(displacements, dtype=np.int32) + 2
- tag = (rank+1)*17 + (i+1)*13
+ tag = (rank+1)*17 + (i+1)*13
tag += dirweight[:directions.size].dot(directions)
tag += disweight[:displacements.size].dot(displacements)
return tag
@@ -302,7 +376,7 @@ def test_mpi_ghost_accumulate(comm):
base_value = ghost_base(rank,directions,displacements,i)
vals = np.full(shape=shape, fill_value=base_value, dtype=dtype)
return vals
-
+
for F in (F0,F1):
dF = F.discretize(topo)
proc_ranks = topo.proc_ranks
@@ -329,12 +403,12 @@ def test_mpi_ghost_accumulate(comm):
for directions in all_directions:
if rank==0:
if directions:
- print ''.join(DirectionLabels[dim-1-d]
+ print ''.join(DirectionLabels[dim-1-d]
for d in directions),
else:
print '--',
sys.stdout.flush()
-
+
for (i,data) in enumerate(dF.data):
data[...] = (rank+1)*(i+1)
for displacements in all_displacements:
@@ -342,16 +416,16 @@ def test_mpi_ghost_accumulate(comm):
continue
(iview,ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
(oview,oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
- data[oview] = ghost_vals(oshape, dtype, rank, directions,
+ data[oview] = ghost_vals(oshape, dtype, rank, directions,
displacements, i)
-
- dF.exchange_ghosts(directions=directions,
+
+ dF.exchange_ghosts(directions=directions,
exchange_method=exchange_method,
ghost_op=GhostOperation.ACCUMULATE)
-
+
for (i,data) in enumerate(dF.data):
for displacements in all_displacements:
- ndisplacements = sum(d!=0 for d in displacements)
+ ndisplacements = sum(d!=0 for d in displacements)
(iview,ishape) = all_inner_ghost_slices[ndirections][directions][displacements]
(oview,oshape) = all_outer_ghost_slices[ndirections][directions][displacements]
@@ -391,12 +465,12 @@ if __name__ == '__main__':
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
-
+
with test_context():
- if (size==0):
+ if (size==1):
test_serial_initialization_1d()
test_serial_initialization_2d()
test_serial_initialization_3d()
-
+
test_mpi_ghost_exchange(comm=comm)
test_mpi_ghost_accumulate(comm=comm)
diff --git a/hysop/mesh/cartesian_mesh.py b/hysop/mesh/cartesian_mesh.py
index 638100860f49ee82a3146e5a8043cd46641bd0e3..2fc95059ffc010ea5873df43e08ad5f749aafe9c 100644
--- a/hysop/mesh/cartesian_mesh.py
+++ b/hysop/mesh/cartesian_mesh.py
@@ -22,11 +22,12 @@ class CartesianMeshView(MeshView):
"""Mesh views on cartesian meshes. """
__slots__ = ('_mesh', '_topology_state')
-
+ isExchInclDiag = True # Is ghosts exchange includes diagonal ghosts
+
@debug
def __new__(cls, mesh, topology_state, **kwds):
"""Create a CartesianMeshView."""
- return super(CartesianMeshView, cls).__new__(cls,
+ return super(CartesianMeshView, cls).__new__(cls,
mesh=mesh, topology_state=topology_state, **kwds)
@debug
@@ -35,22 +36,22 @@ class CartesianMeshView(MeshView):
from hysop.topology.cartesian_topology import CartesianTopologyState
check_instance(mesh, CartesianMesh)
check_instance(topology_state, CartesianTopologyState)
- super(CartesianMeshView, self).__init__(mesh=mesh,
+ super(CartesianMeshView, self).__init__(mesh=mesh,
topology_state=topology_state, **kwds)
def __get_transposed_mesh_attr(self, *attr_names, **kwds):
"""
Get one or more transposed mesh attribute.
- If more than one attribute is requested, return transposed
+ If more than one attribute is requested, return transposed
attributes as a tuple.
- If 'iterable' is set to True (trough kwds), the unique attribute is fetched and
+ If 'iterable' is set to True (trough kwds), the unique attribute is fetched and
its values are transposed afterwards (this implies len(attr_names)==1).
In this case the return value is a tuple.
"""
check_instance(attr_names, tuple, values=str, minsize=1)
-
+
transpose = self._topology_state.transposed
if ('iterable' in kwds) and (kwds['iterable'] is True):
assert len(attr_names)==1
@@ -70,7 +71,7 @@ class CartesianMeshView(MeshView):
def _get_on_proc(self):
"""Return True if the current process has points on the current mesh."""
return self._mesh._on_proc
-
+
def _get_grid_resolution(self):
"""
Effective resolution of the global mesh.
@@ -134,7 +135,7 @@ class CartesianMeshView(MeshView):
def _get_global_boundaries(self):
"""Return global domain boundaries as a tuple of left and right boundaries."""
return self.__get_transposed_mesh_attr('_global_boundaries', iterable=True)
-
+
def _get_local_resolution(self):
"""
Local resolution of this mesh including ghosts.
@@ -217,8 +218,9 @@ class CartesianMeshView(MeshView):
def get_local_inner_ghost_slices(self, ghosts=None):
"""
Return a list of slices defining the ghosts in this arrays as local indices.
- Those slices corresponds to neighbour processes overlap on local
- compute slices (EXCLUDING diagonal processes).
+ Those slices corresponds to neighbour processes overlap on local
+ compute slices (depending on isExchInclDiag parameter it INCLUDES
+ or EXCLUDE diagonal processes).
For each direction, a pair of tuples of slices is returned, one
for left and one for the right ghosts.
"""
@@ -226,20 +228,30 @@ class CartesianMeshView(MeshView):
local_start = self.local_start
local_stop = self.local_stop
compute_resolution = self.compute_resolution
-
+ local_resolution = self.local_resolution
+
ghosts = to_tuple(ghosts or self.ghosts)
ghosts = ghosts*dim if len(ghosts)==1 else ghosts
assert len(ghosts)==dim
+ if self.isExchInclDiag:
+ gh_slices = lambda j: slice(None)
+ else:
+ gh_slices = lambda j: slice(ghosts[j], compute_resolution[j]+ghosts[j])
local_inner_ghost_slices = []
for i in xrange(dim):
- inner_lslices = tuple( slice(ghosts[j], compute_resolution[j]+ghosts[j]) \
- if j!=i else slice(local_start[i],
- local_start[i]+ghosts[i]) for j in xrange(dim) )
- inner_rslices = tuple( slice(ghosts[j], compute_resolution[j]+ghosts[j]) \
- if j!=i else slice(local_stop[i]-ghosts[i],
- local_stop[i]) for j in xrange(dim) )
- inner_shape = compute_resolution.copy()
+ inner_lslices = tuple(
+ gh_slices(j) if j!=i else
+ slice(local_start[i], local_start[i]+ghosts[i])
+ for j in xrange(dim) )
+ inner_rslices = tuple(
+ gh_slices(j) if j!=i else
+ slice(local_stop[i]-ghosts[i], local_stop[i])
+ for j in xrange(dim) )
+ if self.isExchInclDiag:
+ inner_shape = local_resolution.copy()
+ else:
+ inner_shape = compute_resolution.copy()
inner_shape[i] = ghosts[i]
local_inner_ghost_slices.append( (inner_lslices, inner_rslices, inner_shape) )
@@ -248,7 +260,8 @@ class CartesianMeshView(MeshView):
def get_local_outer_ghost_slices(self, ghosts=None):
"""
Return a list of slices defining the ghosts in this arrays as local indices.
- Those slices corresponds to local to process ghosts (EXCLUDING diagonal processes).
+ Those slices corresponds to local to process ghosts (depending
+ on isExchInclDiag it INCLUDES or EXCLUDES diagonal processes).
For each direction, a pair of tuples of slices is returned, one
for left and one for the right ghosts.
"""
@@ -256,34 +269,44 @@ class CartesianMeshView(MeshView):
local_start = self.local_start
local_stop = self.local_stop
compute_resolution = self.compute_resolution
-
+ local_resolution = self.local_resolution
+
ghosts = to_tuple(ghosts or self.ghosts)
ghosts = ghosts*dim if len(ghosts)==1 else ghosts
assert len(ghosts)==dim
+ if self.isExchInclDiag:
+ gh_slices = lambda j: slice(None)
+ else:
+ gh_slices = lambda j: slice(ghosts[j], compute_resolution[j]+ghosts[j])
local_outer_ghost_slices = []
for i in xrange(dim):
- outer_lslices = tuple( slice(ghosts[j], compute_resolution[j]+ghosts[j]) \
- if j!=i else slice(local_start[i]-ghosts[i],
- local_start[i]) for j in xrange(dim) )
- outer_rslices = tuple( slice(ghosts[j], compute_resolution[j]+ghosts[j]) \
- if j!=i else slice(local_stop[i],
- local_stop[i]+ghosts[i]) for j in xrange(dim) )
- outer_shape = compute_resolution.copy()
+ outer_lslices = tuple(
+ gh_slices(j) if j!=i else
+ slice(local_start[i]-ghosts[i], local_start[i])
+ for j in xrange(dim) )
+ outer_rslices = tuple(
+ gh_slices(j) if j!=i else
+ slice(local_stop[i], local_stop[i]+ghosts[i])
+ for j in xrange(dim) )
+ if self.isExchInclDiag:
+ outer_shape = local_resolution.copy()
+ else:
+ outer_shape = compute_resolution.copy()
outer_shape[i] = ghosts[i]
local_outer_ghost_slices.append( (outer_lslices, outer_rslices, outer_shape) )
-
+
return tuple(local_outer_ghost_slices)
-
+
def get_all_local_inner_ghost_slices(self, ghosts=None):
"""
- Return collection of slices and shapes describing all possible
+ Return collection of slices and shapes describing all possible
combinations of inner ghosts slices in this array as local indices.
-
+
Those slices corresponds to local to process ghosts (ie. ghosts that may
be sent to other neighbor processes, including diagonal procecess,
during a ghosts exchange).
-
+
Return a dictionnary {ndirections} (number of directions intersected)
-> {directions} (0=first axis, ..., dim-1=last axis)
-> {displacements} (-1=LEFT, 0=CENTER, +1=RIGHT)
@@ -322,19 +345,19 @@ class CartesianMeshView(MeshView):
shape += (lshape[d]-2*ghosts[d],)
views.setdefault(ndirections, {}) \
.setdefault(directions, {}) \
- .setdefault(displacements, (view,shape))
+ .setdefault(displacements, (view,shape))
return views
-
+
def get_all_local_outer_ghost_slices(self, ghosts=None):
"""
- Return collection of slices and shapes describing all possible combinations
+ Return collection of slices and shapes describing all possible combinations
of outer ghosts slices in this array as local indices.
-
+
Those slices corresponds to local to process ghosts (ie. ghosts that may
be received from other neighbor processes, including diagonal processes,
during a ghosts exchange).
-
+
Return a dictionnary {ndirections} (number of directions intersected)
-> {directions} (0=first axis, ..., dim-1=last axis)
-> {displacements} (-1=LEFT, 0=CENTER, +1=RIGHT)
@@ -373,42 +396,42 @@ class CartesianMeshView(MeshView):
shape += (lshape[d]-2*ghosts[d],)
views.setdefault(ndirections, {}) \
.setdefault(directions, {}) \
- .setdefault(displacements, (view,shape))
+ .setdefault(displacements, (view,shape))
return views
def get_all_local_inner_ghost_slices_per_ncenters(self, ghosts=None):
"""
- Compute the collection of slices describing all possible combinations
- of inner ghosts slice in this array as local indices like
+ Compute the collection of slices describing all possible combinations
+ of inner ghosts slice in this array as local indices like
self.get_all_local_inner_ghost_slices() and sort them by
number of centers (number of displacement == 0).
-
+
Return a list [ndirections] (number of directions intersected)
-> {directions} (0=first axis, ..., dim-1=last axis)
[ncenters] (number of null displacements)
-> {displacements} (-1=LEFT, 0=CENTER, +1=RIGHT)
-> (view,shape)
- -> 'all'
+ -> 'all'
-> tuple of all views and shapes up to this depth
-> ncenters
-> tuple of all views and shapes with given number of centers
"""
views = self.get_all_local_inner_ghost_slices(ghosts=ghosts)
return self._sort_by_centers(views)
-
+
def get_all_local_outer_ghost_slices_per_ncenters(self, ghosts=None):
"""
- Compute the collection of slices describing all possible combinations
- of outer ghosts slice in this array as local indices like
+ Compute the collection of slices describing all possible combinations
+ of outer ghosts slice in this array as local indices like
self.get_all_local_outer_ghost_slices() and sort them by
number of centers (number of displacement == 0).
-
+
Return a dict {ndirections} (number of directions intersected)
-> {directions} (0=first axis, ..., dim-1=last axis)
{ncenters} (number of null displacements)
-> {displacements} (-1=LEFT, 0=CENTER, +1=RIGHT)
-> (view,shape)
- -> 'all'
+ -> 'all'
-> list of all views and shapes up to this depth
-> {ncenters}
-> list of all views and shapes with given number of centers
@@ -441,11 +464,11 @@ class CartesianMeshView(MeshView):
.setdefault(displacements, data)
return views
-
+
def _get_local_origin(self):
"""
- Origin of the first computational point on the physical box
+ Origin of the first computational point on the physical box
including ghosts.
"""
return self.__get_transposed_mesh_attr('_local_origin')
@@ -466,7 +489,7 @@ class CartesianMeshView(MeshView):
Boundaries on the interior of the global domain have value BoundaryCondition.NONE.
"""
return self.__get_transposed_mesh_attr('_local_boundaries', iterable=True)
-
+
def _get_compute_resolution(self):
"""
@@ -493,7 +516,7 @@ class CartesianMeshView(MeshView):
def _get_is_at_right_boundary(self):
"""
Return a numpy boolean mask to identify meshes that are on the right of the domain.
- ie. is_at_right_boundary[d] = True means that process cartesian coordinates
+ ie. is_at_right_boundary[d] = True means that process cartesian coordinates
is the lastest on direction d: proc_coords[d] == proc_shape[d] - 1.
"""
return self.__get_transposed_mesh_attr('_is_at_right_boundary')
@@ -502,14 +525,14 @@ class CartesianMeshView(MeshView):
Return a numpy boolean mask to identify meshes that are on either on the left or on
the right of the domain. Meshes can be on the left and the right at the same time on
direction d if and only if proc_shape[d] == 1.
- ie. is_at_boundary[d] = True means that process cartesian coordinates is the first or
+ ie. is_at_boundary[d] = True means that process cartesian coordinates is the first or
the lastest on direction d:(proc_coords[d] in [0, proc_shape[d] - 1]).
"""
return np.logical_and(self._get_is_at_left_boundary(), self._get_is_at_right_boundary())
-
+
def _get_proc_coords(self):
"""
- Return the coordinate of the process that owns this mesh in the
+ Return the coordinate of the process that owns this mesh in the
CartesianTopology topology.
"""
return self.__get_transposed_mesh_attr('_proc_coords')
@@ -520,7 +543,7 @@ class CartesianMeshView(MeshView):
def iter_axes(self, axes=0):
"""
Return an iterator over a mesh array that returns indices.
- Iterates 'axes' axes at a time (0=whole array view,
+ Iterates 'axes' axes at a time (0=whole array view,
1=one axe at a time, max value is self.dim).
Iterated shape indices are also yielded.
"""
@@ -530,11 +553,11 @@ class CartesianMeshView(MeshView):
for idx in np.ndindex(*iter_shape):
yield (idx, I)
-
+
def iter_compute_mesh(self, axes=0):
"""
Return an iterator over a compute mesh array that returns indices.
- Iterates 'axes' axes at a time (0=whole array view,
+ Iterates 'axes' axes at a time (0=whole array view,
1=one axe at a time, max value is self.dim).
"""
return self.build_compute_mesh_iterator(axes=axes).iter_compute_mesh()
@@ -573,12 +596,12 @@ class CartesianMeshView(MeshView):
def iter_compute_mesh(self):
(dim, ghosts, I, gI) = self._data
- for idx in self._new_ndindex_iterator():
+ for idx in self._new_ndindex_iterator():
gidx = tuple(idx[i]+ghosts[i] for i in xrange(axes))
yield (idx, gidx, I, gI)
return __CartesianMeshComputeAxeIterator(dim=self.dim,
- ghosts=self.ghosts,
+ ghosts=self.ghosts,
compute_resolution=self.compute_resolution,
local_compute_indices=self.local_compute_indices)
@@ -592,7 +615,7 @@ class CartesianMeshView(MeshView):
return False
else:
return True
-
+
def is_inside_global_domain(self, point):
"""
Return true if point is located inside (or on the boundaries)
@@ -606,8 +629,8 @@ class CartesianMeshView(MeshView):
def point_local_indices(self, point):
"""
- Return the local indices of the mesh point that is
- the closest to the given point (which is a position in the
+ Return the local indices of the mesh point that is
+ the closest to the given point (which is a position in the
local subdomain).
If the input point is not inside the local subdomain, return None.
Else return a np.ndarray with dtype HYSOP_DIM.
@@ -618,10 +641,10 @@ class CartesianMeshView(MeshView):
return npw.asintegerarray(indices)
else:
return None
-
+
def point_global_indices(self, point):
"""
- Return the global indices of the mesh point that is
+ Return the global indices of the mesh point that is
the closest to the given point (which is a position in the domain).
If the input point is not inside the domain, return None.
Else return a np.ndarray with dtype HYSOP_DIM.
@@ -632,7 +655,7 @@ class CartesianMeshView(MeshView):
return npw.asintegerarray(indices)
else:
return None
-
+
def local_to_global(self, local_slices):
"""
Convert input local index slices to global index slices.
@@ -652,23 +675,23 @@ class CartesianMeshView(MeshView):
def global_to_local(self, global_slices):
"""
- Perform the intersection of input global mesh indices
- with the local mesh indices of thuis mesh and return
- the resulting global indices that intersected as
+ Perform the intersection of input global mesh indices
+ with the local mesh indices of thuis mesh and return
+ the resulting global indices that intersected as
local indices.
The input is a list or tuple of slices (one per domain direction)
of global indices and the result is a tuple of slices
- of intersected local indices, if the intersection is
+ of intersected local indices, if the intersection is
not empty.
If there is no intersection with the local mesh, return None.
"""
check_instance(global_slices, (list,tuple), values=slice)
-
+
slc = Utils.intersect_slices(global_slices, self.global_compute_slices)
if (slc is None):
- return None
+ return None
else:
return tuple( slice( slc[i].start - self.global_start[i] + self.ghosts[i],
slc[i].stop - self.global_start[i] + self.ghosts[i],
@@ -681,10 +704,10 @@ class CartesianMeshView(MeshView):
"""
s = 'CartesianMesh'
s += '[topo_id={}, proc_coords={}, global_start={}, local_res={}, ghosts={}]'
- s = s.format(self.topology.id, self.proc_coords,
+ s = s.format(self.topology.id, self.proc_coords,
self.global_start, self.local_resolution, self.ghosts)
return s
-
+
def long_description(self):
"""
Long description of this mesh as a string.
@@ -698,7 +721,7 @@ class CartesianMeshView(MeshView):
s += ' *local boundaries: left => {}\n'.format(self.local_boundaries[0])
s += ' right => {}\n'.format(self.local_boundaries[1])
return s
-
+
def __eq__(self, other):
if not isinstance(other, CartesianMeshView):
return NotImplemented
@@ -724,7 +747,7 @@ class CartesianMeshView(MeshView):
topology = property(_get_topology)
on_proc = property(_get_on_proc)
-
+
grid_resolution=property(_get_grid_resolution)
grid_npoints=property(_get_grid_npoints)
@@ -737,7 +760,7 @@ class CartesianMeshView(MeshView):
global_end=property(_get_global_end)
global_length=property(_get_global_length)
global_boundaries=property(_get_global_boundaries)
-
+
local_resolution=property(_get_local_resolution)
local_start=property(_get_local_start)
local_stop=property(_get_local_stop)
@@ -773,14 +796,14 @@ class CartesianMesh(CartesianMeshView, Mesh):
@debug
def __new__(cls, topology, local_resolution, global_start, **kwds):
- return super(CartesianMesh, cls).__new__(cls, topology=topology,
+ return super(CartesianMesh, cls).__new__(cls, topology=topology,
topology_state=topology.topology_state, mesh=None, **kwds)
-
+
@debug
def __init__(self, topology, local_resolution, global_start, **kwds):
"""
A local cartesian grid, defined on each mpi process.
-
+
Parameters
-----------
topology: :class:`hysop.topology.CartesianTopology`
@@ -788,7 +811,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
local_resolution: array like of ints
Local resolution of this mesh, *including* ghost points.
global_start: array like of ints
- Indices in the global mesh of the lowest point of
+ Indices in the global mesh of the lowest point of
the local mesh (excluding ghosts).
kwds: dict
Base class arguments.
@@ -801,7 +824,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
-----
This is a class mainly for internal use not supposed to be called
directly by user but from domain during topology creation.
-
+
Example
-------
Consider a 1D domain and a global resolution with N+1 points
@@ -825,13 +848,13 @@ class CartesianMesh(CartesianMeshView, Mesh):
and you'll get all the local meshes.
- For example, with 2 procs on a 1D periodic domain
+ For example, with 2 procs on a 1D periodic domain
discretized with 9 points:
- global grid (node number): 0 1 2 3 4 5 6 7
+ global grid (node number): 0 1 2 3 4 5 6 7
| | | | | | | |
proc 0 (global indices): X X 0 1 2 3 X X | |
- (local indices) : 0 1 2 3 4 5 6 7 | |
+ (local indices) : 0 1 2 3 4 5 6 7 | |
proc 1 (global indices): X X 4 5 6 7 X X
(local indices): 0 1 2 3 4 5 6 7
with 'X' for ghost points.
@@ -863,17 +886,17 @@ class CartesianMesh(CartesianMeshView, Mesh):
minval=1)
check_instance(global_start, (list,tuple,np.ndarray), values=(int,long,np.integer),
minval=0)
-
- super(CartesianMesh, self).__init__(topology=topology,
+
+ super(CartesianMesh, self).__init__(topology=topology,
topology_state=topology.topology_state, mesh=self, **kwds)
-
+
# ---
# Attributes relative to the global mesh
# i.e common to all MPI processes on the current topology.
# ---
domain = topology.domain
dim = domain.dim
-
+
discretization = topology._topology._discretization
assert (discretization.resolution>=1).all()
assert (discretization.ghosts>=0).all()
@@ -881,12 +904,12 @@ class CartesianMesh(CartesianMeshView, Mesh):
resolution = np.asintegerarray(discretization.resolution).copy()
space_step = npw.asrealarray(domain.length / (resolution - 1))
del resolution, discretization
-
+
global_resolution = npw.asintegerarray(topology.global_resolution).copy()
global_start = npw.asintegerarray(global_start).copy()
assert global_start.size == dim
assert (global_start>=0).all()
-
+
# Remove 1 point on each periodic axe because of periodicity
grid_resolution = global_resolution - domain.periodicity
@@ -899,7 +922,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
assert (local_resolution > 2*ghosts).all()
assert (local_resolution <= grid_resolution + 2*ghosts).all()
compute_resolution = local_resolution - 2*ghosts
-
+
# deduce all other attributes
global_stop = global_start + compute_resolution
global_compute_slices = tuple( slice(i,j) for (i,j) in zip(global_start, global_stop) )
@@ -907,11 +930,11 @@ class CartesianMesh(CartesianMeshView, Mesh):
for i in xrange(dim):
lslices = tuple( slice(None) if j!=i else slice(global_start[i]-ghosts[i], \
global_start[i]) for j in xrange(dim) )
- rslices = tuple( slice(None) if j!=i else slice(global_stop[i],
+ rslices = tuple( slice(None) if j!=i else slice(global_stop[i],
global_stop[i]+ghosts[i]) for j in xrange(dim) )
global_ghost_slices.append( (lslices,rslices) )
global_ghost_slices = tuple(global_ghost_slices)
-
+
global_origin = domain.origin.copy()
global_end = domain.end.copy()
global_length = domain.length.copy()
@@ -919,11 +942,11 @@ class CartesianMesh(CartesianMeshView, Mesh):
local_origin = domain.origin + space_step*(global_start - ghosts)
local_end = domain.origin + space_step*(global_stop + ghosts)
local_length = local_end - local_origin
-
+
local_start = ghosts
local_stop = local_resolution - ghosts
local_compute_slices = tuple( slice(i,j) for (i,j) in zip(local_start, local_stop) )
-
+
local_indices = tuple(np.arange(Ni, dtype=HYSOP_INTEGER) for Ni in local_resolution)
local_coords = tuple(
npw.asrealarray(tuple(local_origin[d] + i*space_step[d]
@@ -931,32 +954,32 @@ class CartesianMesh(CartesianMeshView, Mesh):
for d in xrange(dim))
local_compute_indices = tuple(local_indices[d][local_start[d]:local_stop[d]]
for d in range(dim))
- local_compute_coords = tuple(local_coords[d][local_start[d]:local_stop[d]]
+ local_compute_coords = tuple(local_coords[d][local_start[d]:local_stop[d]]
for d in range(dim))
-
+
# Is this mesh on the last process in some direction in the
# mpi grid of process?
proc_coords = topology.proc_coords
proc_shape = topology.proc_shape
is_at_left_boundary = (proc_coords == 0)
is_at_right_boundary = (proc_coords == proc_shape-1)
-
+
# global boundaries and local boundaries
global_boundaries = domain.boundaries
- lboundaries = np.asarray( [bc if at_left else BoundaryCondition.NONE
+ lboundaries = np.asarray( [bc if at_left else BoundaryCondition.NONE
for (bc, at_left) in zip(global_boundaries[0], is_at_left_boundary) ])
- rboundaries = np.asarray( [bc if at_right else BoundaryCondition.NONE
+ rboundaries = np.asarray( [bc if at_right else BoundaryCondition.NONE
for (bc, at_right) in zip(global_boundaries[1], is_at_right_boundary) ])
local_boundaries = (lboundaries, rboundaries)
-
- npw.set_readonly(global_resolution, grid_resolution, global_start, global_stop,
+
+ npw.set_readonly(global_resolution, grid_resolution, global_start, global_stop,
global_origin, global_end, global_length,
global_boundaries[0], global_boundaries[1],
- local_resolution, local_start, local_stop,
+ local_resolution, local_start, local_stop,
local_origin, local_end, local_length,
- local_boundaries[0], local_boundaries[1],
- compute_resolution, ghosts,
+ local_boundaries[0], local_boundaries[1],
+ compute_resolution, ghosts,
proc_coords, proc_shape,
is_at_left_boundary, is_at_right_boundary,
space_step)
@@ -975,7 +998,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
self._global_end = global_end
self._global_length = global_length
self._global_boundaries = global_boundaries
-
+
self._local_resolution = local_resolution
self._local_start = local_start
self._local_stop = local_stop
@@ -987,9 +1010,9 @@ class CartesianMesh(CartesianMeshView, Mesh):
self._local_coords = local_coords
self._local_compute_indices = local_compute_indices
self._local_compute_coords = local_compute_coords
-
+
self._local_compute_slices = local_compute_slices
-
+
self._compute_resolution = compute_resolution
self._ghosts = ghosts
self._proc_coords = proc_coords
@@ -1021,7 +1044,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
check_instance(self.global_boundaries, tuple, values=np.ndarray, size=2)
for i in xrange(2):
check_instance(self.global_boundaries[i], np.ndarray, dtype=object, size=dim)
-
+
check_instance(self.local_resolution, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
check_instance(self.local_start, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
check_instance(self.local_stop, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
@@ -1052,21 +1075,21 @@ class CartesianMesh(CartesianMeshView, Mesh):
check_instance(self.local_compute_coords, tuple, size=dim, values=np.ndarray)
check_instance(self.local_compute_mesh_coords, tuple, size=dim, values=np.ndarray)
for i in xrange(dim):
- check_instance(self.local_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
+ check_instance(self.local_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
ndim=1, size=self.local_resolution[i])
- check_instance(self.local_mesh_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
+ check_instance(self.local_mesh_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
ndim=dim, size=self.local_resolution[i])
- check_instance(self.local_compute_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
+ check_instance(self.local_compute_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
ndim=1, size=self.compute_resolution[i])
- check_instance(self.local_compute_mesh_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
+ check_instance(self.local_compute_mesh_indices[i], np.ndarray, dtype=HYSOP_INTEGER,
ndim=dim, size=self.compute_resolution[i])
- check_instance(self.local_coords[i], np.ndarray, dtype=HYSOP_REAL,
+ check_instance(self.local_coords[i], np.ndarray, dtype=HYSOP_REAL,
ndim=1, size=self.local_resolution[i])
- check_instance(self.local_mesh_coords[i], np.ndarray, dtype=HYSOP_REAL,
+ check_instance(self.local_mesh_coords[i], np.ndarray, dtype=HYSOP_REAL,
ndim=dim, size=self.local_resolution[dim-1-i])
- check_instance(self.local_compute_coords[i], np.ndarray, dtype=HYSOP_REAL,
+ check_instance(self.local_compute_coords[i], np.ndarray, dtype=HYSOP_REAL,
ndim=1, size=self.compute_resolution[i])
- check_instance(self.local_compute_mesh_coords[i], np.ndarray, dtype=HYSOP_REAL,
+ check_instance(self.local_compute_mesh_coords[i], np.ndarray, dtype=HYSOP_REAL,
ndim=dim, size=self.compute_resolution[dim-1-i])
check_instance(self.compute_resolution, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
@@ -1085,4 +1108,3 @@ class CartesianMesh(CartesianMeshView, Mesh):
from hysop.topology.cartesian_topology import CartesianTopologyState
check_instance(topology_state, CartesianTopologyState)
return CartesianMeshView(mesh=self, topology_state=topology_state)
-