diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index bc49433c9dfdea6ac31c4b4a9781501b0655a4c4..286a08c6d25b9ec334d4693f277c99e8aa11add8 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -799,13 +799,13 @@ class Bridge_intercomm(object):
other_slices[other_rk, 0::2])
intersect_slice[1::2] = np.minimum(my_slices[1::2],
other_slices[other_rk, 1::2])
- # Convert intersection to local slices
+ # Convert intersection to local slices, exluding ghosts
if not (intersect_slice[0::2] > intersect_slice[1::2]).any():
# compatible intersection (else nothing to transfer)
self.transfers[other_rk] = [
- (i - s, j + 1 - s) for i, j, s in zip(
+ (g + i - s, g + j + 1 - s) for i, j, s, g in zip(
intersect_slice[0::2],
- intersect_slice[1::2], my_start)]
+ intersect_slice[1::2], my_start, self.topo.ghosts)]
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 021e9ab10c734903dc22c4eceb2e728c83d1715e..d8d239a140a3c0f6c15a540d4ab5f755c0123c52 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -10,7 +10,7 @@ from parmepy.numerics.differential_operations import Laplacian,\
DivStressTensor3D
from parmepy.numerics.finite_differences import FD_C_2
from parmepy.operator.monitors.monitoring import Monitoring
-from parmepy.mpi import main_comm, main_rank, MPI
+from parmepy.mpi import MPI
from parmepy.tools.timers import timed_function
import parmepy.tools.numpywrappers as npw
import scitools.filetable as ft
@@ -326,7 +326,7 @@ class Forces(Monitoring):
Compute the forces according the Noca s formula
"""
- def __init__(self, velocity, vorticity, topo,
+ def __init__(self, velocity, vorticity, topo,
obstacle, boxMin= None, boxMax=None, Reynolds=None,
uinf=1.0, method=FD_C_2, frequency=None, prefix=None):
"""
@@ -362,7 +362,7 @@ class Forces(Monitoring):
self.output = []
self.topo = self._predefinedTopo
- # Compute coef used to determine drag coefficient :
+ # Compute coef used to determine drag coefficient :
# cD = coef.Fx = 2.Fx/rho.uinf**2.D
self.bufferForce = 0.
rho = 1.0
@@ -385,14 +385,14 @@ class Forces(Monitoring):
self._synchronize = UpdateGhosts(self.topo,
self.velocity.nbComponents)
- # Function for the computation of div(T) where
+ # Function for the computation of div(T) where
# T is the stress tensor: T = Nabla(U) + Nabla(U)t
self.divT = DivStressTensor3D(self.topo, method=self.method)
self.fd_scheme = FD_C_2(self.topo.mesh.space_step)
self._isUpToDate = True
- if (main_rank == 0):
+ if (self.topo.comm.Get_rank() == 0):
self.f = open(self.prefix, 'w')
def compute_control_box(self, obst) :
@@ -470,14 +470,14 @@ class Forces(Monitoring):
# there is no part of the volume control inside the local domain
self.isForceComputationNeeded = False
-# print 'coord_start, coord_end', main_rank, self.coord_start, self.coord_end
-# print 'distMin, distMax, step, res', main_rank, distMin, distMax, self.step, self.res
-# print 'ind_bound_ChiX-', main_rank, ind_boundaries[0][0], coord_boundaries[0][0], self.coords[0][ind_boundaries[0][0], :, :]
-# print 'ind_bound_ChiX+', main_rank, ind_boundaries[0][1], coord_boundaries[0][1], self.coords[0][ind_boundaries[0][1], :, :]
-# print 'ind_bound_ChiY-', main_rank, ind_boundaries[1][0], coord_boundaries[1][0], self.coords[1][:, ind_boundaries[1][0], :]
-# print 'ind_bound_ChiY+', main_rank, ind_boundaries[1][1], coord_boundaries[1][1], self.coords[1][:, ind_boundaries[1][1], :]
-# print 'ind_bound_ChiZ-', main_rank, ind_boundaries[2][0], coord_boundaries[2][0], self.coords[2][:, :, ind_boundaries[2][0]]
-# print 'ind_bound_ChiZ+', main_rank, ind_boundaries[2][1], coord_boundaries[2][1], self.coords[2][:, :, ind_boundaries[2][1]]
+# print 'coord_start, coord_end', self.topo.comm.Get_rank(), self.coord_start, self.coord_end
+# print 'distMin, distMax, step, res', self.topo.comm.Get_rank(), distMin, distMax, self.step, self.res
+# print 'ind_bound_ChiX-', self.topo.comm.Get_rank(), ind_boundaries[0][0], coord_boundaries[0][0], self.coords[0][ind_boundaries[0][0], :, :]
+# print 'ind_bound_ChiX+', self.topo.comm.Get_rank(), ind_boundaries[0][1], coord_boundaries[0][1], self.coords[0][ind_boundaries[0][1], :, :]
+# print 'ind_bound_ChiY-', self.topo.comm.Get_rank(), ind_boundaries[1][0], coord_boundaries[1][0], self.coords[1][:, ind_boundaries[1][0], :]
+# print 'ind_bound_ChiY+', self.topo.comm.Get_rank(), ind_boundaries[1][1], coord_boundaries[1][1], self.coords[1][:, ind_boundaries[1][1], :]
+# print 'ind_bound_ChiZ-', self.topo.comm.Get_rank(), ind_boundaries[2][0], coord_boundaries[2][0], self.coords[2][:, :, ind_boundaries[2][0]]
+# print 'ind_bound_ChiZ+', self.topo.comm.Get_rank(), ind_boundaries[2][1], coord_boundaries[2][1], self.coords[2][:, :, ind_boundaries[2][1]]
# definition of the indicator functions of box boundaries
self.chi_x_minus = [slice(ind_boundaries[0][0], ind_boundaries[0][0] +1),
@@ -570,9 +570,9 @@ class Forces(Monitoring):
localForce = self.integrateOnBox(localForce, vort, self.chi_box, dt)
# Reduction operation
- if main_rank == 0:
+ if self.topo.comm.Get_rank() == 0:
comm = self.topo.comm
- comm.Reduce([localForce, MPI.DOUBLE], [force, MPI.DOUBLE],
+ comm.Reduce([localForce, MPI.DOUBLE], [force, MPI.DOUBLE],
op=MPI.SUM, root=0)
force = force * self.coef
# Print drag, lift and side force coefficients in output file
@@ -616,10 +616,10 @@ class Forces(Monitoring):
## WARNING : this routine is available inly if velo and vort
## topologies have the same global mesh resolution !!!
"""
- Compute integrals on surface normal to "direc" axis to calculate forces
+ Compute integrals on surface normal to "direc" axis to calculate forces
acting on the body. See (2.1) of Noca 1999 or (52) of Plouhmans 2002.
Returns \f$ \int_{S}\mathbf{n}\cdot\gamma\, dS \f$, where
- \f$ \mathbf{n}\cdot\gamma = [(\frac{1}{2} (\mathbf{u}\cdot\mathbf{u})\mathbf{n} -
+ \f$ \mathbf{n}\cdot\gamma = [(\frac{1}{2} (\mathbf{u}\cdot\mathbf{u})\mathbf{n} -
(\mathbf{n}\cdot\mathbf{u})\mathbf{u}) -
\frac{1}{N-1}(\mathbf{n}\cdot\mathbf{u})(\mathbf{x}\wedge\mathbf{\omega}) +
\frac{1}{N-1}(\mathbf{n}\cdot\mathbf{\omega})(\mathbf{x}\wedge\mathbf{u})] +
@@ -648,8 +648,8 @@ class Forces(Monitoring):
self.n_T = self.part3_Zdir
# ====================== part 1 =======================
- # part 1 = 1/2(velocity.velocity)n -
- # (n.velocity)velocity -
+ # part 1 = 1/2(velocity.velocity)n -
+ # (n.velocity)velocity -
# 1/(dim-1)(n.velocity)(coord X vorticity) +
# 1/(dim-1)(n.vorticity)(coord X velocity)
@@ -814,8 +814,3 @@ class Forces(Monitoring):
def __str__(self):
s = "Compute_forces. "
return s
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : Compute_forces"
- print Compute_forces.__doc__
diff --git a/HySoP/hysop/operator/redistribute.py b/HySoP/hysop/operator/redistribute.py
index 211dfd4363883cd1d3bd684751e6e6a6bbc27f18..52ed863cd7a32017fe235db216637e1b4ecfc968 100644
--- a/HySoP/hysop/operator/redistribute.py
+++ b/HySoP/hysop/operator/redistribute.py
@@ -20,7 +20,7 @@ a variable defined on several meshes.
"""
from parmepy import __VERBOSE__
-from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, np
+from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, np, S_DIR
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Bridge
from parmepy.mpi import main_rank, main_comm
@@ -80,6 +80,13 @@ class Redistribute(Operator):
else:
# Only the given component is considered
self._range_components = lambda v: [component]
+ self.r_request = {}
+ self.s_request = {}
+ self._r_types = {}
+ self._s_types = {}
+ for v in self.variables:
+ self._r_types[v] = {}
+ self._s_types[v] = {}
def discretize(self):
pass
@@ -173,6 +180,37 @@ class Redistribute(Operator):
len(self.bridges[v].recvFrom.keys()) == 0 and \
len(self.bridges[v].sendTo.keys()) == 0
+ # Build MPI subarrays
+ dim = self.domain.dimension
+ for v in self.variables:
+ br = self.bridges[v]
+ vTo = self.opTo.discreteFields[v].data[0]
+ vFrom = self.opFrom.discreteFields[v].data[0]
+ for rk in br.recvFrom.keys():
+ subvshape = tuple([br.recvFrom[rk][i].stop -
+ br.recvFrom[rk][i].start
+ for i in range(dim)])
+ substart = tuple([br.recvFrom[rk][i].start
+ for i in range(dim)])
+ self._r_types[v][rk] = \
+ PARMES_MPI_REAL.Create_subarray(vTo.shape,
+ subvshape,
+ substart,
+ order=ORDERMPI)
+ self._r_types[v][rk].Commit()
+ for rk in br.sendTo.keys():
+ subvshape = tuple([br.sendTo[rk][i].stop -
+ br.sendTo[rk][i].start
+ for i in range(dim)])
+ substart = tuple([br.sendTo[rk][i].start
+ for i in range(dim)])
+ self._s_types[v][rk] = \
+ PARMES_MPI_REAL.Create_subarray(vFrom.shape,
+ subvshape,
+ substart,
+ order=ORDERMPI)
+ self._s_types[v][rk].Commit()
+
self._isUpToDate = True
def _apply_toHost_bridges_toDevice(self, simulation=None):
@@ -229,8 +267,7 @@ class Redistribute(Operator):
# (use buffers for mpi send/recv? )
# - move MPI datatypes into the bridge? --> and free MPI type properly
if __VERBOSE__:
- print "{"+str(main_rank)+"}", "_apply_bridges"
- dim = self.domain.dimension
+ print "{" + str(main_rank) + "}", "_apply_bridges"
self.r_request = {}
self.s_request = {}
for v in self.variables:
@@ -238,48 +275,24 @@ class Redistribute(Operator):
# Apply for each component considered
for d in self._range_components(v):
if __VERBOSE__:
- print "{"+str(main_rank)+"}", "_apply_bridges", \
+ print "{" + str(main_rank) + "}", "_apply_bridges", \
self.opFrom.discreteFields[v].name, '->', \
self.opTo.discreteFields[v].name
vTo = self.opTo.discreteFields[v].data[d]
vFrom = self.opFrom.discreteFields[v].data[d]
- v_name = self.opTo.discreteFields[v].name
+ v_name = self.opFrom.discreteFields[v].name + S_DIR[d]
if br.hasLocalInter:
vTo[br.ito] = vFrom[br.ifrom]
for rk in br.recvFrom.keys():
- ind = v_name + str(rk) + '-' + str(d)
- subvshape = tuple([br.recvFrom[rk][i].stop -
- br.recvFrom[rk][i].start
- for i in range(dim)])
- substart = tuple([br.recvFrom[rk][i].start
- for i in range(dim)])
- recvTYPE = PARMES_MPI_REAL.Create_subarray(vTo.shape,
- subvshape,
- substart,
- order=ORDERMPI)
- recvTYPE.Commit()
- self.r_request[ind] = main_comm.Irecv([vTo, 1,
- recvTYPE],
- source=rk,
- tag=rk)
+ self.r_request[v_name + str(rk)] = \
+ main_comm.Irecv([vTo, 1, self._r_types[v][rk]],
+ source=rk, tag=rk)
self._hasRequests = True
for rk in br.sendTo.keys():
- ind = v_name + str(rk) + '-' + str(d)
- subvshape = tuple([br.sendTo[rk][i].stop -
- br.sendTo[rk][i].start
- for i in range(dim)])
- substart = tuple([br.sendTo[rk][i].start
- for i in range(dim)])
- sendTYPE = PARMES_MPI_REAL.Create_subarray(vFrom.shape,
- subvshape,
- substart,
- order=ORDERMPI)
- sendTYPE.Commit()
- self.s_request[ind] = main_comm.Issend([vFrom, 1,
- sendTYPE],
- dest=rk,
- tag=main_rank)
+ self.s_request[v_name + str(rk)] = \
+ main_comm.Issend([vFrom, 1, self._s_types[v][rk]],
+ dest=rk, tag=main_rank)
self._hasRequests = True
def _wait_host(self):
@@ -288,7 +301,7 @@ class Redistribute(Operator):
of all communication requests.
"""
if __VERBOSE__:
- print "{"+str(main_rank)+"}", "WAITING MPI", self._hasRequests
+ print "{" + str(main_rank) + "}", "WAITING MPI", self._hasRequests
if self._hasRequests:
for rk in self.r_request.keys():
self.r_request[rk].Wait()
@@ -298,7 +311,7 @@ class Redistribute(Operator):
def _wait_device(self):
if __VERBOSE__:
- print "{"+str(main_rank)+"}", "WAITING OPENCL"
+ print "{" + str(main_rank) + "}", "WAITING OPENCL"
for dv in self._toDevice_fields + self._toHost_fields:
dv.wait()
@@ -313,11 +326,11 @@ class Redistribute(Operator):
else check for sending to rsend or
receiving from rrecv. Process ranks
should be given in main_comm.
- @param rsend : rank of the process
+ @param rsend : discrete variable name + S_DIR + rank of the process
to which a message has been sent
and for which we want to test
message completion.
- @param rrecv : rank of the process
+ @param rrecv : discrete variable name + S_DIR + rank of the process
from which a message has been receive
and for which we want to test
message completion.
diff --git a/HySoP/hysop/operator/redistribute_intercomm.py b/HySoP/hysop/operator/redistribute_intercomm.py
index 9da815298b883ced57298e6337a056d98ade3b49..a4119ccd8bdc1d201b502a2db11890655feb0226 100644
--- a/HySoP/hysop/operator/redistribute_intercomm.py
+++ b/HySoP/hysop/operator/redistribute_intercomm.py
@@ -7,7 +7,7 @@ the destination.
It relies on a Bridge_intercomm.
"""
-from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI
+from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, S_DIR
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Bridge_intercomm
@@ -58,9 +58,18 @@ class RedistributeIntercomm(Operator):
# Only the given component is considered
self._range_components = lambda v: [component]
+ self._parent_rank = self.parent_comm.Get_rank()
+ self._my_rank = self.topo.comm.Get_rank()
+ self._dim = self.topo.domain.dimension
+
self.bridges = {}
self.r_request = {}
self.s_request = {}
+ self._r_types = {}
+ self._s_types = {}
+ for v in self.variables:
+ self._r_types[v] = {}
+ self._s_types[v] = {}
def discretize(self):
pass
@@ -78,54 +87,102 @@ class RedistributeIntercomm(Operator):
self.id_from, self.id_to,
self.proc_tasks)
+ for v in self.variables:
+ dv = v.discreteFields[self.topo]
+ transfers = self.bridge.transfers
+ # Set reception
+ if self.proc_tasks[self._parent_rank] == self.id_to:
+ for from_rk in transfers.keys():
+ subshape = tuple(
+ [transfers[from_rk][i][1] - transfers[from_rk][i][0]
+ for i in range(self._dim)])
+ substart = tuple(
+ [transfers[from_rk][i][0] for i in range(self._dim)])
+ self._r_types[v][from_rk] = \
+ PARMES_MPI_REAL.Create_subarray(dv[0].shape,
+ subshape,
+ substart,
+ order=ORDERMPI)
+ self._r_types[v][from_rk].Commit()
+ # Set Sending
+ if self.proc_tasks[self._parent_rank] == self.id_from:
+ for to_rk in transfers.keys():
+ subshape = tuple(
+ [transfers[to_rk][i][1] - transfers[to_rk][i][0]
+ for i in range(self._dim)])
+ substart = tuple(
+ [transfers[to_rk][i][0] for i in range(self._dim)])
+ self._r_types[v][to_rk] = \
+ PARMES_MPI_REAL.Create_subarray(dv[0].shape,
+ subshape,
+ substart,
+ order=ORDERMPI)
+ self._r_types[v][to_rk].Commit()
+
def apply(self, simulation=None):
"""
Proceed with data redistribution from opFrom to opTo
"""
- self.r_request = {}
- self.s_request = {}
- parent_rank = self.parent_comm.Get_rank()
- my_rank = self.topo.comm.Get_rank()
- dim = self.topo.domain.dimension
self.parent_comm.Barrier()
+ self.r_request = {}
+ self.s_request = {}
for v in self.variables:
dv = v.discreteFields[self.topo]
transfers = self.bridge.transfers
for d in self._range_components(v):
+ v_name = dv.name + S_DIR[d]
# Set reception
- if self.proc_tasks[parent_rank] == self.id_to:
+ if self.proc_tasks[self._parent_rank] == self.id_to:
for from_rk in transfers.keys():
- subshape = tuple([
- transfers[from_rk][i][1] -
- transfers[from_rk][i][0] for i in range(dim)])
- substart = tuple([
- transfers[from_rk][i][0] for i in range(dim)])
- recv_type = PARMES_MPI_REAL.Create_subarray(
- dv[d].shape, subshape, substart, order=ORDERMPI)
- recv_type.Commit()
- ind = str(from_rk) + '->' + str(my_rank)
- self.r_request[ind] = self.bridge.inter_comm.Irecv(
- [dv[d], 1, recv_type], source=from_rk, tag=from_rk)
+ self.r_request[v_name + str(from_rk)] = \
+ self.bridge.inter_comm.Irecv(
+ [dv[d], 1, self._r_types[v][from_rk]],
+ source=from_rk, tag=from_rk)
# Set Sending
- if self.proc_tasks[parent_rank] == self.id_from:
+ if self.proc_tasks[self._parent_rank] == self.id_from:
for to_rk in transfers.keys():
- subshape = tuple([
- transfers[to_rk][i][1] -
- transfers[to_rk][i][0] for i in range(dim)])
- substart = tuple([
- transfers[to_rk][i][0] for i in range(dim)])
- send_type = PARMES_MPI_REAL.Create_subarray(
- dv[d].shape, subshape, substart, order=ORDERMPI)
- send_type.Commit()
- ind = str(my_rank) + '->' + str(to_rk)
- self.s_request[ind] = self.bridge.inter_comm.Isend(
- [dv[d], 1, send_type], dest=to_rk, tag=my_rank)
+ self.s_request[v_name + str(to_rk)] = \
+ self.bridge.inter_comm.Issend(
+ [dv[d], 1, self._r_types[v][to_rk]],
+ dest=to_rk, tag=self._my_rank)
def wait(self, simulation=None):
"""Wait for requests completion."""
- for rk in self.r_request.keys():
+ for rk in self.r_request:
self.r_request[rk].Wait()
- for rk in self.s_request.keys():
+ for rk in self.s_request:
self.s_request[rk].Wait()
self.parent_comm.Barrier()
+ self.r_request = []
+ self.s_request = []
+
+ def test(self, rsend=None, rrecv=None):
+ """
+ if neither rsend or rrecv is given return
+ True if all communication request are complete
+ else check for sending to rsend or
+ receiving from rrecv. Process ranks
+ should be given in local communicator.
+ @param rsend : variable name + S_DIR + rank
+ @param rrecv : variable name + S_DIR + rank
+ """
+ if(rsend is not None or rrecv is not None):
+ send_res = True
+ recv_res = True
+ if rsend is not None:
+ send_res = self.s_request[rsend].Test()
+ if rrecv is not None:
+ recv_res = self.r_request[rrecv].Test()
+ res = send_res and recv_res
+ else:
+ res = True
+ for rk in self.r_request.keys():
+ res = self.r_request[rk].Test()
+ if not res:
+ return res
+ for rk in self.s_request.keys():
+ res = self.s_request[rk].Test()
+ if not res:
+ return res
+ return res
diff --git a/HySoP/hysop/operator/reprojection.py b/HySoP/hysop/operator/reprojection.py
index c6d0f53552b7cfd183e40fed63e33439976f7cfe..51f6085a3731dc9d68d97d94c3900a356fd86f6d 100644
--- a/HySoP/hysop/operator/reprojection.py
+++ b/HySoP/hysop/operator/reprojection.py
@@ -11,7 +11,7 @@ from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
from parmepy.numerics.differential_operations import GradV
import parmepy.tools.numpywrappers as npw
from parmepy.numerics.updateGhosts import UpdateGhosts
-from parmepy.mpi import main_rank, MPI
+from parmepy.mpi import main_rank, main_comm, main_size, MPI
from parmepy.tools.timers import timed_function