From b85607ba61527a98efb7d4c417c387a138ccd7b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 26 Sep 2014 18:25:55 +0200
Subject: [PATCH] Review mesh mechanism, add/fix subsets. All tests ok except
those with obstacles somewhere. Work in progress
---
HySoP/hysop/domain/mesh.py | 309 +++++++++++
HySoP/hysop/domain/obstacle/controlBox.py | 6 +-
HySoP/hysop/domain/subsets/__init__.py | 3 +
HySoP/hysop/domain/subsets/boxes.py | 156 ++++++
HySoP/hysop/domain/subsets/subset.py | 47 ++
HySoP/hysop/domain/tests/test_mesh.py | 141 +++++
HySoP/hysop/domain/tests/test_submesh.py | 53 ++
HySoP/hysop/domain/tests/test_subset.py | 175 ++++++
HySoP/hysop/mpi/bridge_inter.py | 9 +-
HySoP/hysop/mpi/bridge_overlap.py | 9 +-
HySoP/hysop/mpi/tests/test_mesh.py | 59 --
HySoP/hysop/operator/advold.py | 521 ------------------
HySoP/hysop/operator/computational.py | 2 +-
.../hysop/operator/discrete/drag_and_lift.py | 6 +-
HySoP/hysop/operator/hdf_io.py | 28 +-
HySoP/hysop/tools/misc.py | 44 ++
HySoP/setup.py.in | 1 +
17 files changed, 953 insertions(+), 616 deletions(-)
create mode 100644 HySoP/hysop/domain/mesh.py
create mode 100644 HySoP/hysop/domain/subsets/__init__.py
create mode 100644 HySoP/hysop/domain/subsets/boxes.py
create mode 100644 HySoP/hysop/domain/subsets/subset.py
create mode 100644 HySoP/hysop/domain/tests/test_mesh.py
create mode 100644 HySoP/hysop/domain/tests/test_submesh.py
create mode 100644 HySoP/hysop/domain/tests/test_subset.py
delete mode 100644 HySoP/hysop/mpi/tests/test_mesh.py
delete mode 100644 HySoP/hysop/operator/advold.py
create mode 100644 HySoP/hysop/tools/misc.py
diff --git a/HySoP/hysop/domain/mesh.py b/HySoP/hysop/domain/mesh.py
new file mode 100644
index 000000000..3ca7710fa
--- /dev/null
+++ b/HySoP/hysop/domain/mesh.py
@@ -0,0 +1,309 @@
+"""
+@file mesh.py local cartesian grid.
+"""
+from parmepy.constants import debug
+import parmepy.tools.numpywrappers as npw
+from parmepy.tools.parameters import Discretization
+import numpy as np
+from parmepy.tools.misc import utils
+
+
+class Mesh(object):
+ """
+ A cartesian grid, defined on each mpi process, in relation to a
+ global resolution.
+
+ For example, consider a 1D domain and a global resolution with N+1 points
+ and a ghost layer with 2 points:
+ \code
+ N = 9
+ dom = Box(dimension=1)
+ discr = Discretization([N + 1], [2]
+ \endcode
+
+ \code
+ m = Mesh(dom, discr, [start], [resol])
+ \endcode
+ will describe a grid on the current process, starting at point
+ of index start in the global mesh (from discr),
+ and of local resolution equal to resol.
+
+ Usually, thanks to a topology, a Mesh will be defined
+ on each mpi processes.
+ For example, with 2 procs:
+
+ global grid (node number): 0 1 2 3 4 5 6 7 8
+
+ proc 0 (global indices): X X 0 1 2 3 X X
+ (local indices) : 0 1 2 3 4 5 6 7
+ proc 1 : X X 4 5 6 7 X X
+ 0 1 2 3 4 5 6 7
+ with 'X' for ghost points.
+
+ - Node '8' of the global grid is not represented on local mesh,
+ because of periodicity. N8 == N0
+
+ - on proc 1, we have:
+ - local resolution = 8
+ - global_start = 4
+ - 'computation nodes' = 2, 3, 4, 5
+ - 'ghost nodes' = 0, 1, 6, 7
+
+ Remarks:
+ - all 'global' values refer to the discretization parameter.
+ For example 'global start' = 2 means a point of index 2 in the
+ global resolution.
+ - only periodic grid are considered
+ """
+
+ @debug
+ def __new__(cls, *args, **kw):
+ return object.__new__(cls, *args, **kw)
+
+ @debug
+ def __init__(self, parent, discretization, resolution, global_start):
+ """
+ @param parent : the geometry on which this mesh is defined
+ (it must be a rectangle or a rectangular cuboid)
+ @param discretization : parmepy.tools.parameters.Discretization which
+ describes the global discretization of the domain
+ (resolution and ghost layer)
+ @param global_start
+ """
+ # geometry of reference for the mesh
+ self.domain = parent
+ assert isinstance(discretization, Discretization)
+ ## The global discretization of reference for this mesh.
+ ## It defines the resolution of the grid on the whole domain and
+ ## the number of ghost points on each subdomain.
+ self.discretization = discretization
+
+ ## Local resolution of this mesh, INCLUDING ghost points
+ self.resolution = npw.asdimarray(resolution)
+ self.resolution.flags.writeable = False
+
+ # dimension of the mesh
+ self._dim = self.resolution.size
+
+ ## Position of this mesh in the global grid
+ self.position = []
+ ## Mesh step size
+ self.space_step = self.domain.length / (self.discretization.resolution
+ - 1)
+ # position (index) of the first computational (i.e. no ghosts)
+ # point of this mesh in the global grid.
+ start = npw.asdimarray(global_start)
+
+ # last computational point, global index
+ ghosts = self.discretization.ghosts
+ end = start + self.resolution - 2 * ghosts
+ self.position = [slice(start[d], end[d]) for d in xrange(self._dim)]
+ ## Coordinates of the "lowest" point of this mesh (including ghost)
+ self.origin = self.domain.origin.copy()
+ self.origin += self.space_step * (start - ghosts)
+
+ ## Coordinates of the last point of the mesh
+ self.end = self.origin + self.space_step * (self.resolution - 1)
+
+ # shift between local and global index.
+ # Usefull for start() and end() methods.
+ self._shift = ghosts - start
+
+ ## coordinates of the points of the grid
+ self.coords = None
+
+ ## Local position of the computational points
+ self.iCompute = []
+ # position (index) of the first computational (i.e. no ghosts), local.
+ start = self.discretization.ghosts
+ # last computational point, local index
+ end = self.resolution - start
+ self.iCompute = [slice(start[d], end[d])
+ for d in xrange(self._dim)]
+ self.ind4integ = self.iCompute
+
+ self.coords = np.ix_(*(np.linspace(self.origin[d],
+ self.end[d], self.resolution[d])
+ for d in xrange(self._dim)))
+
+ def local_indices(self, point):
+ """
+ returns indices of the point of coordinates (close to) tab = x, y, z
+ If (x, y, z) is not a grid point, it returns the closest grid point.
+ """
+ point = npw.asrealarray(point)
+ if self.is_inside(point):
+ return npw.asdimarray(np.rint((point - self.origin)
+ / self.space_step))
+ else:
+ return False
+
+ def global_indices(self, point):
+ point = npw.asrealarray(point)
+ return npw.asdimarray(np.rint((point - self.domain.origin)
+ / self.space_step))
+
+ def is_inside(self, point):
+ """
+ @param point: coordinates of a point
+ @return : true if the point is inside the local (mpi proc) mesh
+ """
+ point = npw.asrealarray(point)
+ return ((point - self.origin) >= 0.).all() and ((self.end
+ - point) >= 0.).all()
+
+ def __str__(self):
+ """
+ Sub mesh display
+ """
+ s = 'Local mesh: \n'
+ s += ' - resolution : ' + str(self.resolution) + '\n'
+ s += ' - position in global mesh: ' + str(self.position) + '\n'
+ s += ' - global discretization : ' + str(self.discretization) + '\n'
+ s += ' - computation points :' + str(self.iCompute) + '\n'
+ return s
+
+ def convert2local(self, sl):
+ """
+ Return the local mesh indices from
+ their global indices.
+ Ilocal = Iglobal - _global_start + ghost
+ @param sl: the list of slices (one per dir) for global indices
+ something like [(start_dir1, end_dir1), ...]
+ @return the same kind of list but in local coordinates
+ """
+ tmp = utils.intersl(sl, self.position)
+ if tmp is not None:
+ return [slice(tmp[i].start + self._shift[i],
+ tmp[i].stop + self._shift[i])
+ for i in xrange(self._dim)]
+ else:
+ return None
+
+ def convert2global(self, sl):
+ """
+ Return the global mesh indices from
+ their local indices.
+ Iglobal = Ilocal + _global_start - ghost
+ @param loc_index an array of size domain.dim*2
+ with loc_index = [start_dir1, end_dir1, start_dir2, end_dir2, ...]
+ return an array of indices, same 'setup' as loc_index.
+ """
+ return [slice(sl[i].start - self._shift[i],
+ sl[i].stop - self._shift[i])
+ for i in xrange(self._dim)]
+
+ def __eq__(self, other):
+ """
+ = operator for mesh.
+ We consider that 2 submeshes are equal if all the conditions
+ below are fulfilled:
+ - their global resolution is the same
+ - the ghost layer is the same
+ - their local resolution is the same
+ """
+ if self.__class__ != other.__class__:
+ return False
+ return self.discretization == other.discretization and\
+ npw.equal(self.resolution, other.resolution).all()
+
+ def start(self):
+ """
+ @return indices of the lowest point of the local mesh
+ in the global grid
+ """
+ return npw.asdimarray([p for p in
+ (self.position[d].start
+ for d in xrange(self._dim))])
+
+ def stop(self):
+ """
+ @return indices + 1 of the 'highest' point of the local mesh
+ in the global grid.
+ Warning: made to be used in slices :
+ slice(mesh.start()[d], mesh.stop()[d]) represents the
+ points of the current process mesh, indices given
+ in the global grid.
+ """
+ return npw.asdimarray([p for p in
+ (self.position[d].stop
+ for d in xrange(self._dim))])
+
+
+class SubMesh(object):
+ """
+ A subset of a predefined (distributed) mesh.
+ """
+
+ def __init__(self, mesh, gstart, gend):
+ """
+ @param mesh : the parent mesh
+ @param gstart : indices in the global grid of mesh
+ of the lowest point of this submesh
+ @param gend : indices of the 'highest' point of this submesh,
+ in the global grid.
+ Warning : gend/gstart are global values, that do not depend on the
+ mpi distribution of data.
+
+ \todo : a proper scheme to clarify all the notations for meshes
+ (global/local start end and so on).
+ """
+ ## parent mesh
+ self.mesh = mesh
+ ## Index of the lowest point of the global submesh in the global grid
+ ## of its parent
+ self.gstart = npw.asdimarray(gstart)
+ ## Index of the 'highest' point of the global submesh
+ ## in the global grid of its parent
+ self.gend = npw.asdimarray(gend)
+ ## domain dim
+ self._dim = len(self.gstart)
+ ##
+ self.gposition = [slice(gstart[d], gend[d] + 1)
+ for d in xrange(self._dim)]
+ ##
+ self.discretization = Discretization(self.gend - self.gstart + 1,
+ mesh.discretization.ghosts)
+ dimension = len(self.gstart)
+ # Find which part of submesh is on the current process and
+ # find its computational points. Warning:
+ # the indices of computational points must be
+ # relative to the parent mesh local grid!
+ sl = utils.intersl(self.gposition, self.mesh.position)
+ # Bool to check if this process holds the end (in any direction)
+ # of the domain. Useful for proper integration on this subset.
+ is_last = [False, ] * dimension
+ self.on_proc = sl is not None
+ if self.on_proc:
+ self.position = [s for s in sl]
+ self.iCompute = self.mesh.convert2local(self.position)
+ self.resolution = [self.position[d].stop - self.position[d].start
+ for d in xrange(dimension)]
+ is_last = np.asarray([self.gend[d] <= self.position[d].stop - 1
+ for d in xrange(dimension)])
+ stop = npw.asdimarray([self.iCompute[d].stop
+ for d in xrange(dimension)])
+ stop[is_last] -= 1
+ self.ind4integ = [slice(self.iCompute[d].start,
+ max(stop[d], self.iCompute[d].start + 1))
+ for d in xrange(dimension)]
+ else:
+ self.position = None
+ self.iCompute = [slice(0, 0), ] * dimension
+ self.ind4integ = self.iCompute
+ self.resolution = [0, ] * dimension
+
+ self.coords = [None, ] * dimension
+ self.coords[0] = self.mesh.coords[0][self.iCompute[0]]
+ if dimension > 1:
+ self.coords[1] = self.mesh.coords[1][:, self.iCompute[1], :]
+ if dimension > 2:
+ self.coords[2] = self.mesh.coords[2][:, :, self.iCompute[2]]
+ self.coords = tuple(self.coords)
+ self.coords4int = [None, ] * dimension
+ self.coords4int[0] = self.mesh.coords[0][self.ind4integ[0]]
+ if dimension > 1:
+ self.coords4int[1] = self.mesh.coords[1][:, self.ind4integ[1], :]
+ if dimension > 2:
+ self.coords4int[2] = self.mesh.coords[2][:, :, self.ind4integ[2]]
+ self.coords4int = tuple(self.coords4int)
diff --git a/HySoP/hysop/domain/obstacle/controlBox.py b/HySoP/hysop/domain/obstacle/controlBox.py
index e8588fc7f..9c982da56 100644
--- a/HySoP/hysop/domain/obstacle/controlBox.py
+++ b/HySoP/hysop/domain/obstacle/controlBox.py
@@ -55,7 +55,7 @@ class ControlBox(Obstacle):
## of the obstacle. Only for subspaces.
self.gstart = None
- def createVolumeAndSides(self, spaceStep):
+ def _createVolumeAndSides(self, spaceStep):
"""
@param[in] array of size self._dim, space step size in each direction
This value will be used to compute a tolerance and detect
@@ -112,7 +112,7 @@ class ControlBox(Obstacle):
spaceStep = topo.mesh.space_step
# -- Control volume : union of two halfspaces --
if not self._boxCreated:
- self.createVolumeAndSides(spaceStep)
+ self._createVolumeAndSides(spaceStep)
# Discretize all volume and surfaces of
# the box for topo
@@ -157,7 +157,7 @@ class ControlBox(Obstacle):
# slice(start,end) --> end not included, so +1
end += 1
resol = end - lstart + 2 * topo.mesh.discretization.ghosts
- gstart = lstart + topo.mesh.global_start
+ gstart = lstart + topo.mesh.start()
gstart -= topo.mesh.discretization.ghosts
self.mesh[topo] = SubMesh(self.domain,
topo.mesh.discretization,
diff --git a/HySoP/hysop/domain/subsets/__init__.py b/HySoP/hysop/domain/subsets/__init__.py
new file mode 100644
index 000000000..4dce0d91f
--- /dev/null
+++ b/HySoP/hysop/domain/subsets/__init__.py
@@ -0,0 +1,3 @@
+## @package parmepy.domain.subsets
+# Description (geometry) and discretisation (grid and data distribution) of some subsets
+# of a domain
diff --git a/HySoP/hysop/domain/subsets/boxes.py b/HySoP/hysop/domain/subsets/boxes.py
new file mode 100644
index 000000000..06195693f
--- /dev/null
+++ b/HySoP/hysop/domain/subsets/boxes.py
@@ -0,0 +1,156 @@
+"""
+@file boxes.py
+2D rectangle or 3D Rectangular Cuboid subsets
+"""
+from parmepy.domain.subsets.subset import Subset
+from parmepy.domain.mesh import SubMesh
+import parmepy.tools.numpywrappers as npw
+from parmepy import Field
+import numpy as np
+from parmepy.mpi.topology import Cartesian
+
+
+class SubBox(Subset):
+ """
+ A rectangle (in 2 or 3D space), defined by the coordinates of
+ its lowest point, its lenghts and its normal.
+ """
+ def __init__(self, origin, length, **kwds):
+ """
+ """
+ super(SubBox, self).__init__(**kwds)
+ ## Dictionnary of parmepy.domain.mesh.Submesh, keys=topo, values=mesh
+ self.mesh = {}
+ ## coordinates of the lowest point of this subset
+ self.origin = npw.asrealarray(origin)
+ ## length of this subset
+ self.length = npw.asrealarray(length)
+ ## coordinates of the 'highest' point of this subset
+ self.end = self.origin + self.length
+
+ self._s_dir = np.where(self.length != 0)
+ msg = 'Subset error, the origin is outside of the domain.'
+ assert ((self.origin - self._parent.origin) >= 0).all(), msg
+ msg = 'Subset error, the subset is too large for the domain.'
+ assert ((self._parent.end - self.end) >= 0).all(), msg
+ ## dict of coordinates of the lengthes of the subset (values)
+ ## for a given topo (keys). If the origin/length does not
+ ## fit exactly with the mesh, there may be a small shift of
+ ## these values.
+ self.real_length = {}
+ self.real_orig = {}
+
+ def discretize(self, topo):
+ """
+ Compute a local submesh on the subset, for a given topology
+ @param topo: a parmepy topology
+ """
+ assert isinstance(topo, Cartesian)
+ # Find the indices of starting/ending points in the global_end
+ # grid of the mesh of topo.
+ gstart = topo.mesh.global_indices(self.origin)
+ gend = topo.mesh.global_indices(self.origin + self.length)
+ # create a submesh
+ self.mesh[topo] = SubMesh(topo.mesh, gstart, gend)
+ self.real_length[topo] = (gend - gstart) * topo.mesh.space_step
+ self.real_orig[topo] = topo.domain.origin + \
+ gstart * topo.mesh.space_step
+
+ def apply(self, tab):
+ pass
+
+ def integrate_field_allc(self, field, topo, root=None):
+ """
+ @param[in] field : a parmepy continuous field
+ @param[in] topo : a parmepy.mpi.topology.Cartesian topology
+ @param[in] root : root process used for mpi reduction. If None
+ reduction is done on all processes from topo.
+ @return a numpy array, with res[i] = integral of component i
+ of the input field over the current subset.
+ """
+ res = npw.zeros(field.nbComponents)
+ gres = npw.zeros(field.nbComponents)
+ for i in xrange(res.size):
+ res[i] = self._integrate_field_on_proc(field, topo, component=i)
+ if root is None:
+ topo.comm.Allreduce(res, gres)
+ else:
+ topo.comm.Reduce(res, gres, root=root)
+
+ return gres
+
+ def integrate_field(self, field, topo, component=0, root=None):
+ """
+ @param[in] field : a parmepy continuous field
+ @param[in] topo : a parmepy.mpi.topology.Cartesian topology
+ @param[in] component : number of the component of field
+ to integrate
+ @return integral of a component of the input field
+ over the current subset, for the discretization given by
+ topo
+ """
+ res = self._integrate_field_on_proc(field, topo, component)
+ if root is None:
+ return topo.comm.allreduce(res)
+ else:
+ return topo.comm.reduce(res, root=root)
+
+ def _integrate_field_on_proc(self, field, topo, component=0):
+ """
+ @param[in] field : a parmepy continuous field
+ @param[in] topo : a parmepy.mpi.topology.Cartesian topology
+ @param[in] component : number of the component of field
+ integrate the field on the subset on the current process
+ (i.e. no mpi reduce over all processes) for the discretization
+ given by topo.
+ @return integral of a component of the input field
+ over the current subset.
+ """
+ assert isinstance(field, Field)
+ assert isinstance(topo, Cartesian)
+ df = field.discretize(topo)
+ dvol = npw.prod(topo.mesh.space_step[self._s_dir])
+ ic = self.mesh[topo].ind4integ
+ result = npw.real_sum(df[component][ic])
+ result *= dvol
+ return result
+
+ def integrate_func(self, func, topo, nbc, root=None):
+ """
+ @param[in] func : a python function, depending on the space
+ coordinates. See below for the required signature.
+ @param[in] topo : a parmepy.mpi.topology.Cartesian topology
+ @param[in] nbc : number of components of the return value from func
+ @return integral of a component of the input field
+ over the current subset, for the discretization given by
+ topo
+ """
+ res = npw.zeros(nbc)
+ gres = npw.zeros(nbc)
+ for i in xrange(res.size):
+ res[i] = self._integrate_func_on_proc(func, topo)
+ if root is None:
+ topo.comm.Allreduce(res, gres)
+ else:
+ topo.comm.Reduce(res, gres, root=root)
+ return gres
+
+ def _integrate_func_on_proc(self, func, topo):
+ """
+ @param[in] func : a python function, depending on the space
+ coordinates. See below for the required signature.
+ @param[in] : a parmepy.mpi.topology.Cartesian topology
+ integrate a funcion on the subset on the current process
+ (i.e. no mpi reduce over all processes)
+ """
+ assert hasattr(func, '__call__')
+ assert isinstance(topo, Cartesian)
+ coords = self.mesh[topo].coords4int
+ vfunc = np.vectorize(func)
+ if self.mesh[topo].on_proc:
+ result = npw.real_sum(vfunc(*coords))
+ else:
+ result = 0.
+ dvol = npw.prod(topo.mesh.space_step)
+ result *= dvol
+ return result
diff --git a/HySoP/hysop/domain/subsets/subset.py b/HySoP/hysop/domain/subsets/subset.py
new file mode 100644
index 000000000..d9d749768
--- /dev/null
+++ b/HySoP/hysop/domain/subsets/subset.py
@@ -0,0 +1,47 @@
+"""
+@file subset.py
+Abstract interface to subset of a given domain
+"""
+from abc import ABCMeta, abstractmethod
+from parmepy.domain.domain import Domain
+
+
+class Subset(object):
+ """
+ A subset is a geometry defined inside a domain.
+ Given a topology on the parent domain, the subset must
+ provide 'slices' or 'tabs' to allow the computation of a discrete
+ field on the subset, with:
+
+ data[subset.slices] = ...
+ or
+ data[subset.tab] = ...
+ or
+ subset.apply(data) = ...
+
+ slices : a list of python slices defining the grid points inside the subset
+ tab : an array of boolean, where tab[i, j, k] = True
+ if the grid point at (i, j, k) is inside the subset.
+
+ There are two types of subsets:
+ - 'regular' ones, those on which a regular mesh can be defined
+ - others, like spheres, cylinders ...
+
+ For regular subsets, one must use 'slices' to compute field inside
+ the subset whereas 'tab' is required for others.
+ """
+
+ __metaclass__ = ABCMeta
+
+ def __init__(self, parent):
+ """
+ @param parent : the domain in which the subset is defined
+ """
+ assert isinstance(parent, Domain)
+ self._parent = parent
+
+ @abstractmethod
+ def apply(self, field):
+ """
+ return field on the current subset
+ """
diff --git a/HySoP/hysop/domain/tests/test_mesh.py b/HySoP/hysop/domain/tests/test_mesh.py
new file mode 100644
index 000000000..0999bd537
--- /dev/null
+++ b/HySoP/hysop/domain/tests/test_mesh.py
@@ -0,0 +1,141 @@
+"""
+@file parmepy.domain.tests.test_mesh
+Testing regular grids.
+"""
+from parmepy.domain.box import Box
+from parmepy.tools.parameters import Discretization
+from parmepy.mpi import main_size, main_rank
+import numpy as np
+
+
+Nx = Ny = Nz = 32
+
+
+def create_mesh(discr):
+ """
+ Periodic mesh
+ """
+ gh = discr.ghosts
+ dim = gh.size
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ cutdir = [False, ] * dim
+ cutdir[-1] = True
+ topo = dom.create_topology(discr, cutdir=cutdir)
+ mesh = topo.mesh
+ # Test global grid
+ assert mesh.discretization == discr
+ assert (mesh.space_step == dom.length / (discr.resolution - 1)).all()
+ # Position compared with global grid
+ shift = np.asarray([0, ] * dim)
+ shift[-1] = Nz / main_size * main_rank
+ resolution = (discr.resolution - 1) / topo.shape + 2 * gh
+ end = shift + resolution - 2 * gh
+ assert mesh.position == [slice(shift[d], end[d]) for d in xrange(dim)]
+ assert (mesh.start() == shift).all()
+ assert (mesh.stop() == end).all()
+ # Test local grid
+ assert (mesh.resolution == resolution).all()
+ assert (mesh.origin == dom.origin + (shift - gh) * mesh.space_step).all()
+ assert mesh.iCompute == [slice(gh[d], resolution[d] - gh[d])
+ for d in xrange(dim)]
+ assert (mesh.global_indices(dom.origin) == [0, ] * dim).all()
+ assert (mesh.local_indices(mesh.origin) == [0, ] * dim).all()
+ if main_size == 0:
+ assert mesh.is_inside([0.3, ] * dim)
+ pt2 = [0., ] * dim
+ pt2[-1] = - 0.8
+ assert not mesh.is_inside(pt2)
+ point = 6 * mesh.space_step + dom.origin
+ point[-1] = (Nz - 1) * mesh.space_step[-1] + dom.origin[-1]
+ req_point = [6, ] * dim
+ req_point[-1] = Nz - 1
+ assert (mesh.global_indices(point) == req_point).all()
+ if main_rank == main_size - 1:
+ assert mesh.is_inside(point)
+ pos = np.asarray(req_point)
+ pos[-1] = Nz / main_size - 1
+ pos += gh
+ assert (mesh.local_indices(point) == pos).all()
+ else:
+ assert mesh.local_indices(point) is False
+
+
+def test_mesh3D():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1])
+ create_mesh(discr)
+
+
+def test_mesh3D_ghost():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [2, 2, 2])
+ create_mesh(discr)
+
+
+def test_mesh2D():
+ discr = Discretization([Nx + 1, Nz + 1])
+ create_mesh(discr)
+
+
+def test_mesh2D_ghost():
+ discr = Discretization([Nx + 1, Nz + 1], [2, 2])
+ create_mesh(discr)
+
+
+def test_convert_local():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [2, 1, 2])
+ gh = discr.ghosts
+ dim = gh.size
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ mesh = topo.mesh
+
+ # test conversion
+ start = 2
+ end = 10
+ source = [slice(start, end), ] * dim
+ res = mesh.convert2local(source)
+ if main_size == 1:
+ assert res == [slice(start + gh[d], end + gh[d]) for d in xrange(dim)]
+ elif main_size == 8:
+ newstart = start - topo.rank * Nz / main_size
+ newend = end - topo.rank * Nz / main_size
+ slref = slice(max(mesh.iCompute[-1].start, newstart + gh[-1]),
+ min(max(newend + gh[-1], gh[-1]),
+ mesh.iCompute[-1].stop))
+ sl = [slice(start + gh[d], end + gh[d]) for d in xrange(dim)]
+ sl[-1] = slref
+ if res is not None:
+ assert res == sl
+ else:
+ assert slref.stop - slref.start == 0
+
+
+def test_convert_global():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [2, 1, 2])
+ gh = discr.ghosts
+ dim = gh.size
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ mesh = topo.mesh
+
+ # test conversion
+ start = 3
+ end = 6
+ source = [slice(start, end), ] * dim
+ res = mesh.convert2global(source)
+ if main_size == 1:
+ assert res == [slice(start - gh[d], end - gh[d]) for d in xrange(dim)]
+ elif main_size == 8:
+ newstart = start + topo.rank * Nz / main_size
+ newend = end + topo.rank * Nz / main_size
+ slref = slice(newstart - gh[-1], newend - gh[-1])
+ sl = [slice(start - gh[d], end - gh[d]) for d in xrange(dim)]
+ sl[-1] = slref
+ assert res == sl
+
+if __name__ == '__main__':
+ test_mesh3D()
+ test_mesh3D_ghost()
+ test_mesh2D()
+ test_mesh2D_ghost()
+ test_convert_local()
+ test_convert_global()
diff --git a/HySoP/hysop/domain/tests/test_submesh.py b/HySoP/hysop/domain/tests/test_submesh.py
new file mode 100644
index 000000000..f3135b6eb
--- /dev/null
+++ b/HySoP/hysop/domain/tests/test_submesh.py
@@ -0,0 +1,53 @@
+"""
+@file parmepy.domain.tests.test_mesh
+Testing regular grids.
+"""
+from parmepy.domain.box import Box
+from parmepy.tools.parameters import Discretization
+from parmepy.domain.mesh import SubMesh
+from parmepy.tools.misc import utils
+
+Nx = Ny = Nz = 32
+g = 2
+
+
+def test_submesh():
+ """
+ Periodic mesh
+ """
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ gh = discr.ghosts
+ dim = gh.size
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ mesh = topo.mesh
+ # Set a starting point and a length for the submesh
+ hh = mesh.space_step
+ xr = 1. * hh
+ # Find the indices of these points
+ gstart = mesh.global_indices(xr)
+ gend = mesh.global_indices(xr + 15 * hh)
+ # Create the submesh
+ subm = SubMesh(mesh, gstart, gend)
+ # check global params of the submesh
+ assert subm.mesh is mesh
+ assert (subm.gstart == gstart).all()
+ assert (subm.gend == gend).all()
+ # the global position we expect
+ gpos = [slice(gstart[d], gend[d] + 1) for d in xrange(dim)]
+ assert subm.gposition == gpos
+ assert (subm.discretization.resolution == [gend[d] - gstart[d] + 1
+ for d in xrange(dim)]).all()
+ assert (subm.discretization.ghosts == mesh.discretization.ghosts).all()
+ # check position of the local submesh, relative to the global grid
+ # It must corresponds to the intersection of the expected position
+ # of the submesh and the position of the parent mesh.
+ if subm.on_proc:
+ assert subm.position == [s for s in utils.intersl(gpos, mesh.position)]
+
+ r2 = mesh.convert2local([s for s
+ in utils.intersl(gpos, mesh.position)])
+ assert subm.iCompute == r2
+
+if __name__ == '__main__':
+ test_submesh()
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
new file mode 100644
index 000000000..e0a291679
--- /dev/null
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -0,0 +1,175 @@
+from parmepy.domain.subsets.boxes import SubBox
+from parmepy.tools.parameters import Discretization
+from parmepy import Field, Box
+import numpy as np
+import math
+
+
+Nx = Ny = Nz = 128
+g = 2
+
+
+def v3d(res, x, y, z, t):
+ res[0][...] = 1.
+ res[1][...] = 1.
+ res[2][...] = 1.
+ return res
+
+
+def f_test(x, y, z):
+ return 1
+
+
+def f_test_2(x, y, z):
+ return math.cos(z)
+
+
+def test_subbox():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ # Starting point and length of the subdomain
+ xr = [0.2, 0.4, 0.1]
+ lr = [0.3, 0.4, 0.8]
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ assert (cub.length == lr).all()
+ assert (cub.origin == xr).all()
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ cub.discretize(topo)
+ assert cub.mesh.values()[0].mesh == topo.mesh
+ assert cub.mesh.keys()[0] == topo
+
+
+def test_integ():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.01, ] * dim)
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ hh = topo.mesh.space_step
+ # Starting point and length of the subdomain
+ xr = 3 * hh
+ lr = 10 * hh
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ cub.discretize(topo)
+ velo = Field(domain=dom, isVector=True, formula=v3d, name='velo')
+ velo.discretize(topo)
+ velo.initialize(topo=topo)
+ i0 = cub.integrate_field_allc(velo, topo)
+ vref = np.prod(cub.real_length[topo])
+ assert (np.abs(i0 - vref) < 1e-6).all()
+
+
+def test_integ_2():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.01, ] * dim)
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ # Starting point and length of the subdomain
+ xr = [0.1, 0.2, 0.15]
+ lr = [0.51, 0.53, 0.3]
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ cub.discretize(topo)
+ velo = Field(domain=dom, isVector=True, formula=v3d, name='velo')
+ velo.discretize(topo)
+ velo.initialize(topo=topo)
+ vref = np.prod(cub.real_length[topo])
+ for i in xrange(velo.nbComponents):
+ i0 = cub.integrate_field(velo, topo, component=i)
+ assert np.abs(i0 - vref) < 1e-6
+
+
+def test_integ_3():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.01, ] * dim)
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ # Starting point and length of the subdomain
+ xr = [0.1, 0.2, 0.15]
+ lr = [0.51, 0.53, 0.3]
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ cub.discretize(topo)
+ nbc = 1
+ vref = np.prod(cub.real_length[topo])
+ for i in xrange(1, nbc + 1):
+ i0 = cub.integrate_func(f_test, topo, nbc=i)
+ assert np.abs(i0 - vref) < 1e-6
+
+
+def test_integ_4():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0., ] * dim,
+ length=[math.pi * 2., ] * dim)
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ # Starting point and length of the subdomain
+ xr = [math.pi * 0.5, ] * dim
+ lr = [math.pi, ] * dim
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ cub.discretize(topo)
+ nbc = 1
+ vref = - 2 * np.prod(cub.real_length[topo][:2])
+ for i in xrange(1, nbc + 1):
+ i0 = cub.integrate_func(f_test_2, topo, nbc=i)
+ assert np.abs(i0 - vref) / np.abs(vref) < topo.mesh.space_step[2]
+
+
+def test_subbox_2d():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ # Starting point and length of the subdomain
+ xr = [0.2, 0.4, 0.1]
+ lr = [0., 0.2, 0.8]
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ assert (cub.length == lr).all()
+ assert (cub.origin == xr).all()
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ cub.discretize(topo)
+ assert cub.mesh.values()[0].mesh == topo.mesh
+ assert cub.mesh.keys()[0] == topo
+ velo = Field(domain=dom, isVector=True, formula=v3d, name='velo')
+ velo.discretize(topo)
+ velo.initialize(topo=topo)
+ i0 = cub.integrate_field_allc(velo, topo)
+ vref = np.prod(cub.real_length[topo][1:])
+ assert (np.abs(i0 - vref) < 1e-6).all()
+
+
+def test_line():
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ dim = 3
+ dom = Box(dimension=dim, origin=[0.1, ] * dim)
+ # Starting point and length of the subdomain
+ xr = [0.2, 0.4, 0.1]
+ lr = [0., 0., 0.8]
+ cub = SubBox(origin=xr, length=lr, parent=dom)
+ assert (cub.length == lr).all()
+ assert (cub.origin == xr).all()
+ # discretization of the dom and of its subset
+ topo = dom.create_topology(discr, cutdir=[False, False, True])
+ cub.discretize(topo)
+ assert cub.mesh.values()[0].mesh == topo.mesh
+ assert cub.mesh.keys()[0] == topo
+ velo = Field(domain=dom, isVector=True, formula=v3d, name='velo')
+ velo.discretize(topo)
+ velo.initialize(topo=topo)
+ i0 = cub.integrate_field_allc(velo, topo)
+ vref = cub.real_length[topo][2]
+ assert (np.abs(i0 - vref) < 1e-6).all()
+
+
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+ test_subbox()
+ test_integ()
+ test_integ_2()
+ test_integ_3()
+ test_subbox_2d()
+ test_line()
diff --git a/HySoP/hysop/mpi/bridge_inter.py b/HySoP/hysop/mpi/bridge_inter.py
index ae500e18f..55fd153a7 100644
--- a/HySoP/hysop/mpi/bridge_inter.py
+++ b/HySoP/hysop/mpi/bridge_inter.py
@@ -4,6 +4,7 @@ Tools to compute the intersection between
two parmes topologies.
"""
from parmepy.mpi.topology import Cartesian, topotools
+from parmepy.tools.misc import utils
from parmepy.mpi import MPI
import parmepy.tools.numpywrappers as npw
@@ -72,12 +73,12 @@ class BridgeInter(object):
self._tranfer_indices = {}
current = current_indices[self._topology.rank]
for rk in remote_indices:
- inter = topotools.intersl(current, remote_indices[rk])
+ inter = utils.intersl(current, remote_indices[rk])
if inter is not None:
self._tranfer_indices[rk] = inter
# Back to local indices
- convert = self._topology.mesh.toIndexLocal2
+ convert = self._topology.mesh.convert2local
self._tranfer_indices = {rk: convert(self._tranfer_indices[rk])
for rk in self._tranfer_indices}
@@ -115,8 +116,8 @@ class BridgeInter(object):
self.comm.bcast(current_indices, root=MPI.PROC_NULL)
# Convert numpy arrays to dict of slices ...
- current_indices = topotools.arrayToDict(current_indices)
- remote_indices = topotools.arrayToDict(remote_indices)
+ current_indices = utils.arrayToDict(current_indices)
+ remote_indices = utils.arrayToDict(remote_indices)
return current_indices, remote_indices
diff --git a/HySoP/hysop/mpi/bridge_overlap.py b/HySoP/hysop/mpi/bridge_overlap.py
index d03f87911..cdc1007a1 100644
--- a/HySoP/hysop/mpi/bridge_overlap.py
+++ b/HySoP/hysop/mpi/bridge_overlap.py
@@ -5,6 +5,7 @@ two parmes topologies defined on the same comm but for a
different number of processes.
"""
from parmepy.mpi.topology import Cartesian, topotools
+from parmepy.tools.misc import utils
from parmepy.mpi.bridge import Bridge
@@ -89,13 +90,13 @@ class BridgeOverlap(Bridge):
current = [slice(None, None, None), ] * dimension
for rk in indices_target:
- inter = topotools.intersl(current, indices_target[rk])
+ inter = utils.intersl(current, indices_target[rk])
if inter is not None:
self._send_indices[rk] = inter
if self._source is not None:
# Back to local indices
- convert = self._source.mesh.toIndexLocal2
+ convert = self._source.mesh.convert2local
self._send_indices = {rk: convert(self._send_indices[rk])
for rk in self._send_indices}
@@ -104,11 +105,11 @@ class BridgeOverlap(Bridge):
else:
current = [slice(None, None, None), ] * dimension
for rk in indices_source:
- inter = topotools.intersl(current, indices_source[rk])
+ inter = utils.intersl(current, indices_source[rk])
if inter is not None:
self._recv_indices[rk] = inter
if self._target is not None:
- convert = self._target.mesh.toIndexLocal2
+ convert = self._target.mesh.convert2local
self._recv_indices = {rk: convert(self._recv_indices[rk])
for rk in self._recv_indices}
diff --git a/HySoP/hysop/mpi/tests/test_mesh.py b/HySoP/hysop/mpi/tests/test_mesh.py
deleted file mode 100644
index ac416debf..000000000
--- a/HySoP/hysop/mpi/tests/test_mesh.py
+++ /dev/null
@@ -1,59 +0,0 @@
-"""
-@file parmepy.mpi.tests.test_mesh
-Testing mesh.
-"""
-from parmepy.domain.box import Box
-from parmepy.tools.parameters import Discretization
-
-
-def test_mesh3D():
- """Periodic mesh"""
- dom = Box()
- resolTopo = Discretization([33, 33, 17])
- topo = dom.create_topology(resolTopo, dim=3)
- res = resolTopo.resolution
- dx = [1. / (n - 1) for n in res]
- assert (topo.mesh.origin == [0. for n in res]).all()
- assert (topo.mesh.end == [1. - dxx for dxx in dx]).all()
- assert (topo.mesh.global_start == [0 for n in res]).all()
- assert (topo.mesh.global_end == [n - 2 for n in res]).all()
- assert (topo.mesh.local_start == [0 for n in res]).all()
- assert (topo.mesh.local_end == [n - 2 for n in res]).all()
- assert len(topo.mesh.coords) == 3
- assert topo.mesh.coords[0].shape == (res[0] - 1, 1, 1)
- assert topo.mesh.coords[1].shape == (1, res[1] - 1, 1)
- assert topo.mesh.coords[2].shape == (1, 1, res[2] - 1)
- assert topo.mesh.coords[0][1, 0, 0] == dx[0]
- assert topo.mesh.coords[1][0, 1, 0] == dx[1]
- assert topo.mesh.coords[2][0, 0, 1] == dx[2]
-
-
-def test_mesh3D_ghost():
- """Periodic mesh"""
- dom = Box()
- resolTopo = Discretization([33, 33, 17], [2, 2, 2])
- topo = dom.create_topology(resolTopo, dim=3)
- ghost = topo.ghosts()
- res = resolTopo.resolution
- dx = [1. / (n - 1) for n in res]
- assert (topo.mesh.origin ==
- [0. - g * dxx for dxx, g in zip(dx, ghost)]).all()
- assert (topo.mesh.end ==
- [1. - dxx + g * dxx for dxx, g in zip(dx, ghost)]).all()
- assert (topo.mesh.global_start == [0 for n in res]).all()
- assert (topo.mesh.global_end == [n - 2 for n in res]).all()
- assert (topo.mesh.local_start == [g for g in ghost]).all()
- assert (topo.mesh.local_end ==
- [n - 2 + g for n, g in zip(res, ghost)]).all()
- assert len(topo.mesh.coords) == 3
- assert topo.mesh.coords[0].shape == (res[0] - 1 + 2 * ghost[0], 1, 1)
- assert topo.mesh.coords[1].shape == (1, res[1] - 1 + 2 * ghost[1], 1)
- assert topo.mesh.coords[2].shape == (1, 1, res[2] - 1 + 2 * ghost[2])
- assert topo.mesh.coords[0][1, 0, 0] == (-ghost[0] + 1) * dx[0]
- assert topo.mesh.coords[1][0, 1, 0] == (-ghost[1] + 1) * dx[1]
- assert topo.mesh.coords[2][0, 0, 1] == (-ghost[2] + 1) * dx[2]
-
-# Todo : update tests for multi-proc mpi runs.
-#if __name__ == '__main__':
-# test_mesh3D()
-# test_mesh3D_ghost()
diff --git a/HySoP/hysop/operator/advold.py b/HySoP/hysop/operator/advold.py
deleted file mode 100644
index 7980d0e9b..000000000
--- a/HySoP/hysop/operator/advold.py
+++ /dev/null
@@ -1,521 +0,0 @@
-"""
-@file advection.py
-
-Advection of a field.
-"""
-from parmepy.constants import debug, np, PARMES_INDEX, S_DIR
-from parmepy.operator.computational import Computational
-from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
- Remesh, Support, Splitting, MultiScale
-from parmepy.numerics.remeshing import L2_1
-from parmepy.fields.continuous import Field
-from parmepy.operator.redistribute import Redistribute
-from parmepy.operator.advectionDir import AdvectionDir
-import parmepy.default_methods as default
-import parmepy.tools.numpywrappers as npw
-from parmepy.operator.continuous import opsetup
-
-
-class Advection(Computational):
- """
- Advection of a field,
- \f{eqnarray*}
- X = Op(X,velocity)
- \f} for :
- \f{eqnarray*}
- \frac{\partial{X}}{\partial{t}} + velocity.\nabla X = 0
- \f}
- Note : we assume incompressible flow.
-
- Computations are performed within a dimensional splitting as folows:
- - 2nd order:
- - X-dir, half time step
- - Y-dir, half time step
- - Z-dir, full time step
- - Y-dir, half time step
- - X-dir, half time step
- - 2nd order full half-steps:
- - X-dir, half time step
- - Y-dir, half time step
- - Z-dir, half time step
- - Z-dir, half time step
- - Y-dir, half time step
- - X-dir, half time step
- - 1st order g:
- - X-dir, half time step
- - Y-dir, half time step
- - Z-dir, half time step
-
- """
-
- @debug
- def __init__(self, velocity, advectedFields, **kwds):
- """
- Create a Transport operator from given variables velocity and scalar.
-
- @param velocity : velocity variable.
- @param advectedFields : Advected fields (may be a list of Fields).
- @param resolutions : list of resolutions (one per variable)
- @param method : Method used
- @param splittingConfig : Dimensional splitting configuration
- (default 'o2')
- @param topo : a predefined topology to discretize variables
- """
- v = [velocity]
- if isinstance(advectedFields, list):
- self.advectedFields = advectedFields
- else:
- self.advectedFields = [advectedFields]
- v += self.advectedFields
- assert 'variables' not in kwds, 'variables parameter is useless.'
- super(Advection, self).__init__(variables=v, **kwds)
-
- vars_str = "_("
- for vv in self.advectedFields:
- vars_str += vv.name + ","
- self.name += vars_str[0:-1] + ')'
-
- if self.method is None:
- self.method = default.ADVECTION
-
- ## Transport velocity
- self.velocity = velocity
- ## Transported fields
- self.output = self.advectedFields
- self.input = [var for var in self.variables]
-
- self.config = {}
-
- # Find which solver is used for advection,
- # among Scales, pure-python and GPU-like.
- # Check also operator-splitting type.
- if Scales in self.method.keys():
- self._isSCALES = True
- # Default splitting = Strang
- if Splitting not in self.method.keys():
- self.method[Splitting] = 'strang'
- else:
- self._isSCALES = False
- assert TimeIntegrator in self.method.keys()
- assert Interpolation in self.method.keys()
- assert Remesh in self.method.keys()
- assert Support in self.method.keys()
- if Splitting not in self.method.keys():
- self.method[Splitting] = 'o2'
- self._old_dir = 0
- self.splitting = []
- self._dim = self.velocity.dimension
- self.advecDir = None
- if not self._isSCALES:
- particles_advectedFields = [
- Field(adF.domain, name="Particle_AdvectedFields",
- isVector=adF.isVector)
- for adF in self.advectedFields]
- if self.method[Support].find('gpu_1k') >= 0:
- particles_positions = None
- else:
- particles_positions = \
- Field(self.advectedFields[0].domain,
- name="Particle_Position", isVector=False
- )
-
- # Directional continuous Advection operators
- self.advecDir = [None] * self._dim
- for i in xrange(self._dim):
- self.advecDir[i] = AdvectionDir(
- self.velocity, self.advectedFields, i,
- particles_advectedFields, particles_positions,
- isMultiScale=self._isMultiScale,
- name_suffix=vars_str[0:-1] + ')', **kwds)
-
- # function to switch between CPU or GPU setup.
- self._my_setup = None
-
- @debug
- def discretize(self):
- """
- Discretisation according to the chosen method.
- Available methods :
- - 'scales' : SCALES fortran routines (3d only, list of vector
- and/or scalar)
- - 'p_O2' : order 4 method, corrected to allow large CFL number,
- untagged particles
- - 'p_O4' : order 4 method, corrected to allow large CFL number,
- untagged particles
- - 'p_L2' : limited and corrected lambda 2
- - 'p_M4' : Lambda_2,1 (=M'4) 4 point formula
- - 'p_M6' (default) : Lambda_4,2 (=M'6) 6 point formula
- - 'p_M8' : M8prime formula
- - 'p_44' : Lambda_4,4 formula
- - 'p_64' : Lambda_6,4 formula
- - 'p_66' : Lambda_6,6 formula
- - 'p_84' : Lambda_8,4 formula
- - 'gpu' : OpenCL kernels (2d and 3d, single field, scalar or vector)
- - Kernels versions:
- - '1k' : Single OpenCL kernel for advection and remeshing
- - '2k' : Separate kernels
- - Integration method:
- - 'rk2' : Runge Kutta 2nd order advection
- - 'rk4' : Runge Kutta 4th order advection
- - remeshing formula:
- - 'm4prime' : = 'l2_1'
- - 'l2_1' : Labmda2,1 : (=M'4) 4 point formula, C1 regularity
- - 'l2_2' : Labmda2,2 : 4 point formula, C2 regularity
- - 'm6prime' : = 'l4_2'
- - 'l4_2' : Labmda4,2 : (=M'6) 6 point formula, C2 regularity
- - 'l4_3' : Labmda4,3 : 6 point formula, C3 regularity
- - 'l4_4' : Labmda4,4 : 6 point formula, C4 regularity
- - 'l6_3' : Labmda6,3 : 8 point formula, C3 regularity
- - 'l6_4' : Labmda6,4 : 8 point formula, C4 regularity
- - 'l6_5' : Labmda6,5 : 8 point formula, C5 regularity
- - 'l6_6' : Labmda6,6 : 8 point formula, C6 regularity
- - 'l8_4' : Labmda8,4 : 10 point formula, C4 regularity
- - 'm8prime' : M8prime formula
- - other : Pure python (2d and 3d, list of vector and/or scalar)
- - Integration method:
- - 'rk2' : Runge Kutta 2nd order advection
- - 'rk4' : Runge Kutta 4th order advection
- - remeshing formula:
- - 'm4prime' : = 'l2_1'
- - 'l2_1' : Labmda2,1 : (=M'4) 4 point formula, C1 regularity
- - 'l2_2' : Labmda2,2 : 4 point formula, C2 regularity
- - 'm6prime' : = 'l4_2'
- - 'l4_2' : Labmda4,2 : (=M'6) 6 point formula, C2 regularity
- - 'l4_3' : Labmda4,3 : 6 point formula, C3 regularity
- - 'l4_4' : Labmda4,4 : 6 point formula, C4 regularity
- - 'l6_3' : Labmda6,3 : 8 point formula, C3 regularity
- - 'l6_4' : Labmda6,4 : 8 point formula, C4 regularity
- - 'l6_5' : Labmda6,5 : 8 point formula, C5 regularity
- - 'l6_6' : Labmda6,6 : 8 point formula, C6 regularity
- - 'l8_4' : Labmda8,4 : 10 point formula, C4 regularity
- - 'm8prime' : M8prime formula
- """
- # --- Advection solver from SCALES ---
- if self._isSCALES:
- if not self._dim == 3:
- raise ValueError("Scales Advection not implemented in 2D.")
- # - Scales imports -
- from parmepy.f2py import scales2py as scales
-
- # - Extract order form self.method (default p_M6) -
- order = None
- for o in ['p_O2', 'p_O4', 'p_L2',
- 'p_M4', 'p_M6', 'p_M8',
- 'p_44', 'p_64', 'p_66', 'p_84']:
- if self.method[Scales].find(o) >= 0:
- order = o
- if order is None:
- print ('Unknown advection method, turn to default (p_M6).')
- order = 'p_M6'
- # - Extract splitting form self.method (default strang) -
- splitting = 'strang'
- for s in ['classic', 'strang', 'particle']:
- if self.method[Splitting].find(s) >= 0:
- splitting = s
-
- # - Create the topologies (get param from scales) -
- # Scales nbcells equals resolutions - 1
- nbcells = npw.asintarray(self.resolutions[self.advectedFields[0]])
- nbcells -= 1
- if self._fromTopo:
- main_size = self.topology.size
- comm = self.topology.comm
- topodims = self.topology.shape
- elif self._comm is not None:
- main_size = self._comm.Get_size()
- comm = self._comm
- topodims = [1, 1, main_size]
- else:
- from parmepy.mpi.main_var import main_size
- from parmepy.mpi.main_var import main_comm as comm
- topodims = [1, 1, main_size]
- scalesres, scalesoffset = \
- scales.init_advection_solver(nbcells,
- self.domain.length,
- topodims, comm.py2f(),
- order=order,
- dim_split=splitting)
- msg = 'Scales Advection not yet implemented with ghosts points.'
- # Use same topodims as scales to create Cartesian topology
- # in order to discretize our fields
- # Case 1 : topology provided by user at init
- # and one topo for all fields
- if self._fromTopo:
- for v in self.variables:
- topo = self.topologies[v]
- assert (topo.shape == topodims).all(),\
- 'input topology is not scales compliant.'
- assert not (topo.ghosts > 0).any(), msg
- self.discreteFields[v] = v.discretize(topo)
-
- # Case 2 : a dict of resolutions
- else:
- if self.ghosts is not None:
- assert not (self.ghosts > 0).any(), msg
- if self._singleTopo:
- topo = self.domain.getOrCreateTopology(
- self._dim, self.resolutions[self.velocity],
- topodims, precomputed=True, offset=scalesoffset,
- localres=scalesres, ghosts=self.ghosts,
- comm=self._comm)
- for v in self.variables:
- self.discreteFields[v] = v.discretize(topo)
- else:
- for v in self.variables:
- topo = \
- self.domain.getOrCreateTopology(
- self._dim, self.resolutions[v], topodims,
- precomputed=True, offset=scalesoffset,
- localres=scalesres, ghosts=self.ghosts,
- comm=self._comm)
- # ... and the corresponding discrete field
- self.discreteFields[v] = v.discretize(topo)
-
-
- if self._isMultiScale:
- self.config['isMultiscale'] = self._isMultiScale
- v_shape = np.asarray(self.resolutions[self.velocity],
- dtype=PARMES_INDEX) - 1
- scales.init_multiscale(v_shape[0], v_shape[1], v_shape[2],
- self.method[MultiScale])
- self._my_setup = self.setup_Scales
-
- # --- GPU or pure-python advection ---
- else:
- for i in xrange(self._dim):
- self.advecDir[i].discretize()
- self.discreteFields = self.advecDir[0].discreteFields
- self._my_setup = self.setup_Python
-
- @staticmethod
- def getWorkLengths(method=None, domain_dim=None):
- """
- Return the length of working arrays lists required
- for advction discrete operator, depending on :
- - the time integrator (RK2, ...)
- - the interpolation (which depends on domain dimension)
- - the remeshing (which depends on domain dimension)
- @param method : the dict of parameters for the operator.
- Default = parmepy.default_methods.ADVECTION
- """
- if method is None:
- method = default.ADVECTION
- assert Interpolation in method,\
- 'An interpolation is required for the advection method.'
- assert TimeIntegrator in method,\
- 'A time integrator is required for the advection method.'
- assert Remesh in method,\
- 'A remesh is required for the advection method.'
- tw = method[TimeIntegrator].getWorkLengths(1)
- iw, iiw = method[Interpolation].getWorkLengths(domain_dim=domain_dim)
- rw, riw = method[Remesh].getWorkLengths(domain_dim=domain_dim)
- return max(tw + iw, rw), max(iiw, riw)
-
- def setWorks(self, rwork=None, iwork=None):
- if rwork is None:
- rwork = []
- if iwork is None:
- iwork = []
- if not self._isSCALES:
- for i in xrange(self._dim):
- self.advecDir[i].setWorks(rwork, iwork)
-
- @opsetup
- def setup(self):
- # Check resolutions to set multiscale case, if required.
- self._isMultiScale = False
- v_resol = self.variables[self.velocity]
- if v_resol != self.variables[self.advectedFields[0]]:
- self._isMultiScale = True
- if self._isMultiScale and not MultiScale in self.method.keys():
- print ("Using default mutiscale interpolation : L2_1")
- self.method[MultiScale] = L2_1
-
- if not self._is_uptodate:
- self._my_setup()
-
- def setup_Scales(self):
- advectedDiscreteFields = [self.discreteFields[f]
- for f in self.advectedFields]
- # - Create the discreteOperator from the
- # list of discrete fields -
- from parmepy.operator.discrete.scales_advection import \
- ScalesAdvection
- self.discreteOperator = ScalesAdvection(
- self.discreteFields[self.velocity],
- advectedDiscreteFields, method=self.method,
- **self.config)
-
- # -- Final set up --
- self.discreteOperator.setup()
- self._is_uptodate = True
-
- def setup_Python(self):
- for i in xrange(self._dim):
- self.advecDir[i].setup()
- # If topologies differs between directions,
- # one need Redistribute
- # operators
- main_size = self.advecDir[0].discreteFields[
- self.velocity].topology.size
- if main_size > 1:
- # Build bridges
- self.bridges = self._dim * [None]
- if main_size > 1:
- for dfrom in xrange(self.domain.dimension):
- self.bridges[dfrom] = self._dim * [None]
- for dto in xrange(self._dim):
- if dfrom == dto:
- self.bridges[dfrom][dto] = None
- else:
- nsuffix = str(dfrom) + '_' + str(dto)
- self.bridges[dfrom][dto] = Redistribute(
- variables=self.advectedFields,
- opFrom=self.advecDir[dfrom],
- opTo=self.advecDir[dto],
- name_suffix=nsuffix)
- self.bridges[dfrom][dto].setup()
-
- # Splitting configuration
- if self.method[Splitting] == 'o2_FullHalf':
- ## Half timestep in all directions
- [self.splitting.append((i, 0.5))
- for i in xrange(self._dim)]
- [self.splitting.append((self._dim - 1 - i, 0.5))
- for i in xrange(self._dim)]
- elif self.method[Splitting] == 'o1':
- [self.splitting.append((i, 1.)) for i in xrange(self._dim)]
- elif self.method[Splitting] == 'x_only':
- self.splitting.append((0, 1.))
- elif self.method[Splitting] == 'y_only':
- self.splitting.append((1, 1.))
- elif self.method[Splitting] == 'z_only':
- self.splitting.append((2, 1.))
- elif self.method[Splitting] == 'o2':
- ## Half timestep in all directions but last
- [self.splitting.append((i, 0.5))
- for i in xrange(self._dim - 1)]
- self.splitting.append((self._dim - 1, 1.))
- [self.splitting.append((self._dim - 2 - i, 0.5))
- for i in xrange(self._dim - 1)]
- else:
- raise ValueError('Unknown splitting configuration:' +
- self.method[Splitting])
-
- if self.method[Support].find('gpu') >= 0:
- splitting_nbSteps = len(self.splitting)
- for d in xrange(self.domain.dimension):
- dOp = self.advecDir[d].discreteOperator
- assert len(dOp.exec_list) == splitting_nbSteps, \
- "Discrete operator execution " + \
- "list and splitting steps sizes must be equal " + \
- str(len(dOp.exec_list)) + " != " + \
- str(splitting_nbSteps)
-
- if main_size > 1:
- self.apply = self._apply_Comm
- else:
- self.apply = self._apply_noComm
- self._is_uptodate = True
-
- if self.method[Support].find('gpu') >= 0:
- s = ""
- device_id = self.advecDir[d].discreteOperator.cl_env._device_id
- gpu_comm = self.advecDir[d].discreteOperator.cl_env.gpu_comm
- gpu_rank = gpu_comm.Get_rank()
- if gpu_rank == 0:
- s += "=== OpenCL buffers allocated"
- s += " on Device:{0} ===\n".format(device_id)
- s += "Global memory used:\n"
- total_gmem = 0
- for d in xrange(self.domain.dimension):
- g_mem_d = 0
- for df in self.advecDir[d].discreteOperator.variables:
- if not df.gpu_allocated:
- df.allocate()
- g_mem_df = gpu_comm.allreduce(df.mem_size)
- g_mem_d += g_mem_df
- if gpu_rank == 0:
- s += " Advection" + S_DIR[d] + ": {0:9d}".format(g_mem_d)
- s += "Bytes ({0:5d} MB)\n".format(g_mem_d / (1024 ** 2))
- total_gmem += g_mem_d
- if gpu_rank == 0:
- s += " Total : {0:9d}".format(total_gmem)
- s += "Bytes ({0:5d} MB)\n".format(total_gmem / (1024 ** 2))
- s += "Local memory used:\n"
- total_lmem = 0
- for d in xrange(self.domain.dimension):
- l_mem_d = gpu_comm.allreduce(
- self.advecDir[d].discreteOperator.size_local_alloc)
- if gpu_rank == 0:
- s += " Advection" + S_DIR[d] + ": {0:9d}".format(l_mem_d)
- s += "Bytes ({0:5d} MB)\n".format(l_mem_d / (1024 ** 2))
- total_lmem += l_mem_d
- if gpu_rank == 0:
- s += " Total : {0:9d}".format(total_lmem) + "Bytes"
- print (s)
-
- @debug
- def _apply_noComm(self, simulation=None):
- """
- Apply this operator to its variables.
- @param simulation : object that describes the simulation
- parameters (time, time step, iteration number ...), see
- parmepy.problem.simulation.Simulation for details.
-
- Redefinition for advection. Applying a dimensional splitting.
- """
- for req in self.requirements:
- req.wait()
- for split_id, split in enumerate(self.splitting):
- self.advecDir[split[0]].apply(
- simulation, split[1], split_id, self._old_dir)
- self._old_dir = split[0]
-
- @debug
- def _apply_Comm(self, simulation=None):
- """
- Apply this operator to its variables.
- @param simulation : object that describes the simulation
- parameters (time, time step, iteration number ...), see
- parmepy.problem.simulation.Simulation for details.
-
- Redefinition for advection. Applying a dimensional splitting.
- """
- for req in self.requirements:
- req.wait()
- for split_id, split in enumerate(self.splitting):
- # Calling the redistribute operators between directions
- if not self._old_dir == split[0]:
- self.bridges[self._old_dir][split[0]].apply()
- self.bridges[self._old_dir][split[0]].wait()
- self.advecDir[split[0]].apply(
- simulation, split[1], split_id, self._old_dir)
- self._old_dir = split[0]
-
- @debug
- def finalize(self):
- """
- Memory cleaning.
- """
- if self._isSCALES:
- Computational.finalize(self)
- else:
- for dop in self.advecDir:
- dop.finalize()
- self.timer = self.timer + dop.timer
-
- def __str__(self):
- """
- Common printings for operators
- """
- shortName = str(self.__class__).rpartition('.')[-1][0:-2]
- if self._isSCALES:
- super(Advection, self).__str__()
- else:
- for i in xrange(self._dim):
- if self.advecDir[i].discreteOperator is not None:
- s = str(self.advecDir[i].discreteOperator)
- else:
- s = shortName + " operator. Not discretised."
- return s + "\n"
diff --git a/HySoP/hysop/operator/computational.py b/HySoP/hysop/operator/computational.py
index 4a51b4cc8..80b30340c 100644
--- a/HySoP/hysop/operator/computational.py
+++ b/HySoP/hysop/operator/computational.py
@@ -264,7 +264,7 @@ class Computational(Operator):
resolution, self.domain.length, comm=comm.py2f())
assert (topo.mesh.resolution == localres).all()
- assert (topo.mesh.global_start == global_start).all()
+ assert (topo.mesh.start() == global_start).all()
msg = 'Ghosts points not yet implemented for fftw-type operators.'
assert (topo.ghosts() == 0).all(), msg
diff --git a/HySoP/hysop/operator/discrete/drag_and_lift.py b/HySoP/hysop/operator/discrete/drag_and_lift.py
index 0d54f4284..a22a48a75 100644
--- a/HySoP/hysop/operator/discrete/drag_and_lift.py
+++ b/HySoP/hysop/operator/discrete/drag_and_lift.py
@@ -10,7 +10,7 @@ from parmepy.numerics.finite_differences import FD_C_2
from parmepy.operator.discrete.discrete import DiscreteOperator
import parmepy.tools.numpywrappers as npw
from parmepy.tools.profiler import profile
-from parmepy.numerics.utils.Utils import sum_cross_product
+from parmepy.numerics.utils.Utils import Utils
class DragAndLift(DiscreteOperator):
@@ -242,7 +242,7 @@ class DragAndLift(DiscreteOperator):
# slices representing the points inside the subset
# (global indices relative to the global mesh of the current topology)
sl = self._voc.slices[self._topology]
- res = sum_cross_product(coords, self.vorticity.data,
- sl, self._rwork, self._obsinds)
+ res = Utils.sum_cross_product(coords, self.vorticity.data,
+ sl, self._rwork, self._obsinds)
res *= self._dvol
return res
diff --git a/HySoP/hysop/operator/hdf_io.py b/HySoP/hysop/operator/hdf_io.py
index 25715432e..6338a081f 100644
--- a/HySoP/hysop/operator/hdf_io.py
+++ b/HySoP/hysop/operator/hdf_io.py
@@ -37,7 +37,7 @@ class HDF_IO(Computational):
@param var_names : a dictionnary of names to connect fields
to the dataset in the hdf file. See example below.
@param subset : a subset of the domain, on which data are read.
- It must be a parmepy.domain.obstacles.Obstacle.
+ It must be a parmepy.domain.subset.
Names paramater example:
if variables=[velo, vorti], and if hdf file contains
@@ -108,26 +108,13 @@ class HDF_IO(Computational):
# Discretize the subset, if required
if self.subset is not None:
self.subset.discretize(self._topology)
- self._slices = self.subset.slices[self._topology]
- # Global resolution for hdf5 output
- self._global_resolution = \
- self.subset.global_resolution(self._topology)
- g_start = self.subset.gstart
- # convert self._slices to global position in topo
- sl = self._topology.mesh.toIndexGlobal(self._slices)
- # And shift using global position of the surface
- sl = [slice(sl[i].start - g_start[i], sl[i].stop - g_start[i])
- for i in xrange(self.domain.dimension)]
-
+ refmesh = self.subset.mesh[self._topology]
else:
- self._global_resolution = \
- list(self._topology.mesh.discretization.resolution - 1)
- self._slices = self._topology.mesh.iCompute
- g_start = self._topology.mesh.global_start
- g_end = self._topology.mesh.global_end + 1
- sl = [slice(g_start[i], g_end[i])
- for i in xrange(self.domain.dimension)]
-
+ refmesh = self._topology.mesh
+ # Global resolution for hdf5 output
+ self._global_resolution = list(refmesh.discretization.resolution - 1)
+ self._slices = refmesh.iCompute
+ sl = list(refmesh.position)
# Reverse order, to fit with xdmf req.
self._global_resolution.reverse()
sl.reverse()
@@ -136,7 +123,6 @@ class HDF_IO(Computational):
@opsetup
def setup(self, rwork=None, iwork=None):
# No list of hdf dataset names provided by user ...
- print 'HDF setup ...', self.__class__
if self.var_names is None:
# Get field names and initialize dataset dict.
for df in self.discreteFields.values():
diff --git a/HySoP/hysop/tools/misc.py b/HySoP/hysop/tools/misc.py
new file mode 100644
index 000000000..dca1db598
--- /dev/null
+++ b/HySoP/hysop/tools/misc.py
@@ -0,0 +1,44 @@
+"""
+@file misc.py
+Some utilities to deal with slices and other python objects
+"""
+
+
+class utils(object):
+ @staticmethod
+ def arrayToDict(inarray):
+ """
+ convert an array into a dictionnary,
+ keys being the column numbers in array
+ and values the content of each corresponding column
+ transformed into a list of slices like this:
+ column = [1, 4, 2, 6, ...] ---> [slice(1, 4), slice(2, 6), ...]
+ """
+ outslice = {}
+ size = inarray.shape[1]
+ dimension = (int)(0.5 * inarray.shape[0])
+ for rk in xrange(size):
+ outslice[rk] = [slice(inarray[2 * d, rk],
+ inarray[2 * d + 1, rk] + 1)
+ for d in xrange(dimension)]
+ return outslice
+
+ @staticmethod
+ def intersl(sl1, sl2):
+ """
+ @param sl1 : a list of slices
+ @param sl2 : a list of slices
+ @return : guess what ... a list of slices such that:
+ result[i] = intersection between sl1[i] and sl2[i]
+ """
+ assert len(sl1) == len(sl2)
+ res = [None] * len(sl1)
+ for d in xrange(len(sl1)):
+ s1 = sl1[d]
+ s2 = sl2[d]
+ start = max(s1.start, s2.start)
+ stop = min(s1.stop, s2.stop)
+ if stop <= start:
+ return None
+ res[d] = slice(start, stop)
+ return tuple(res)
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 04c5da8e6..86e72c608 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -14,6 +14,7 @@ name = '@PYPACKAGE_NAME@'
# List of modules (directories) to be included
packages = ['parmepy',
'parmepy.domain',
+ 'parmepy.domain.subsets',
'parmepy.domain.obstacle',
'parmepy.fields',
'parmepy.operator',
--
GitLab