Fix last bugs in subsets. All tests ok except python advection and...

Fix last bugs in subsets. All tests ok except python advection and penalization (--> obstacles obsolete)

HySoP/hysop/domain/subsets/control_box.py

+"""
+@file control_box.py
+Define a volume of control inside a domain (volume + all faces)
+"""
+from parmepy.domain.subsets.boxes import SubBox
+import numpy as np
+import parmepy.tools.numpywrappers as npw
+
+
+class ControlBox(SubBox):
+    """
+    Build a sub-domain, box-shaped
+    ==> define set of indices inside this domain (ind member)
+    and set of indices belonging to surfaces of this domain (slices members).
+    Useful to define control volume to perform integration.
+    """
+
+    def __init__(self, **kwds):
+        """
+        Build the volume of control
+        @param origin : coordinates of the lowest point in the sub-domain
+        @param lengths : lengths of box sides.
+        """
+        super(ControlBox, self).__init__(**kwds)
+
+        # We must have a real box, not a plane ...
+        assert len(self.s_dir) == len(self.length)
+
+        self.surf = [None] * len(self.length) * 2
+
+    def discretize(self, topo):
+        # Create the mesh for the whole box
+        super(ControlBox, self).discretize(topo)
+        # Create a mesh for each side
+        dim = topo.domain.dimension
+        ilist = np.arange(dim)
+        for direction in xrange(dim):
+            ndir = np.where(ilist == direction)[0]
+            length = self.real_length[topo].copy()
+            length[ndir] = 0.0
+            orig = self.real_orig[topo].copy()
+            self.surf[2 * direction] = \
+                SubBox(origin=orig, length=length, parent=self._parent)
+            orig = self.real_orig[topo].copy()
+            orig[ndir] += self.real_length[topo][ndir]
+            self.surf[2 * direction + 1] = \
+                SubBox(origin=orig, length=length, parent=self._parent)
+
+        for s in self.surf:
+            s.discretize(topo)
+
+    def _check_boundaries(self, surf, topo):
+        """
+        Special care when some boundaries of the control box are on the
+        upper boundaries of the domain.
+        Remind that for periodic bc, such surfaces does not really
+        exists in the parent mesh.
+        """
+        return surf.mesh[topo].check_boundaries()
+
+    def integrate_on_faces(self, field, topo, list_dir,
+                           component=0, 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.
+        @param[in] list_dir : list of surfaces on which we must integrate.
+        0 : normal to xdir, lower surf,
+        1 : normal to xdir, upper surf (in x dir)
+        2 : normal to ydir, lower surf, and so on.
+        the field must be integrated.
+        @return a numpy array, with res = sum
+        of the integrals of a component of field on all surf in list_dir
+        """
+        res = 0.
+        for ndir in list_dir:
+            surf = self.surf[ndir]
+            msg = 'Control Box error : surface out of bounds.'
+            assert self._check_boundaries(surf, topo), msg
+            res += surf.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_on_faces_allc(self, field, topo, list_dir, 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.
+        @param[in] list_dir : list of surfaces on which we must integrate.
+        0 : normal to xdir, lower surf,
+        1 : normal to xdir, upper surf (in x dir)
+        2 : normal to ydir, lower surf, and so on.
+        the field must be integrated.
+        @return a numpy array, with res[i] = sum
+        of the integrals of component i of field on all surf in list_dir
+        """
+        nbc = field.nbComponents
+        res = npw.zeros(nbc)
+        gres = npw.zeros(nbc)
+        for ndir in list_dir:
+            surf = self.surf[ndir]
+            assert self._check_boundaries(surf, topo)
+            for i in xrange(nbc):
+                res[i] += surf.integrate_field_on_proc(field, topo, i)
+        if root is None:
+            topo.comm.Allreduce(res, gres)
+        else:
+            topo.comm.Reduce(res, gres, root=root)
+        return gres

HySoP/hysop/domain/subsets/submesh.py

+"""
+@file submesh.py local cartesian grid.
+"""
+import parmepy.tools.numpywrappers as npw
+from parmepy.tools.parameters import Discretization
+import numpy as np
+from parmepy.tools.misc import utils
+
+
+class SubMesh(object):
+    """
+    A subset of a predefined (distributed) mesh.
+    """
+
+    def __init__(self, mesh, substart, subend):
+        """
+        @param mesh : the parent mesh
+        @param substart : indices in the global grid of mesh
+        of the lowest point of this submesh
+        @param subend : indices of the 'highest' point of this submesh,
+        in the global grid.
+        Warning : subend/substart 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).
+        """
+        # Note : all variables with 'global' prefix are related
+        # to the global grid, that is the mesh defined on the whole
+        # domain. These variables do not depend on the mpi distribution.
+
+        # ---
+        # Attributes relative to the global mesh
+
+        ## parent mesh
+        self.mesh = mesh
+        # dimension of the submesh
+        self._dim = self.mesh.domain.dimension
+
+        ## Index of the lowest point of the global submesh in the global grid
+        ## of its parent
+        self.substart = npw.asdimarray(substart)
+        ## Index of the 'highest' point of the global submesh
+        ## in the global grid of its parent
+        self.subend = npw.asdimarray(subend)
+
+        # position of the submesh in the global grid of its parent mesh.
+        global_position_in_parent = [slice(substart[d], self.subend[d] + 1)
+                                     for d in xrange(self._dim)]
+        hh = self.mesh.space_step
+        ## Coordinates of the lowest point of this submesh
+        self.global_origin = self.substart * hh + self.mesh.domain.origin
+        ## Length of this submesh
+        self.global_length = (self.subend - self.substart) * hh
+
+        # Warning : we must not overpass the parent global discretization.
+        gres = self.subend - self.substart + 1
+        # directions where length is 0, i.e. directions 'normal' to
+        # the submesh.
+        self._n_dir = np.where(gres == 1)[0]
+        ## discretization of the subset
+        ## Warning : at the time, no ghosts on the submesh!
+        self.discretization = Discretization(gres)
+        # 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(global_position_in_parent, 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.
+        self.is_last = [False, ] * self._dim
+
+        ## Check if a part of the submesh is present on the current proc.
+        self.on_proc = sl is not None
+
+        if self.on_proc:
+            # Is this mesh on the last process in some direction in the
+            # mpi grid of process?
+            self.is_last = np.asarray([self.subend[d] < sl[d].stop
+                                       for d in xrange(self._dim)])
+            ## position of the LOCAL submesh in the global grid
+            ## of the parent mesh
+            self.position_in_parent = [s for s in sl]
+            ## Indices of the points of the submesh, relative to
+            ## the LOCAL array
+            self.iCompute = self.mesh.convert2local(self.position_in_parent)
+
+            ## Resolution of the local submesh
+            self.resolution = [self.iCompute[d].stop - self.iCompute[d].start
+                               for d in xrange(self._dim)]
+
+            # Now, we must find which points must be used
+            # when we integrate on this submesh
+            stops = npw.asdimarray([self.iCompute[d].stop
+                                    for d in xrange(self._dim)])
+            # ilist = [d for d in ll if self.is_last[d]]
+            # stops[ilist] += 1
+            # when 'is_last', the last point must be removed for integration
+            stops[self.is_last] -= 1
+            # and finally, for direction where subset length is zero,
+            # we must increase stop, else integral will always be zero!
+            stops[self._n_dir] = npw.asdimarray([self.iCompute[d].start + 1
+                                                 for d in self._n_dir])
+
+            ## Same as self.iCompute but recomputed to be used
+            ## for integration on the submesh
+            self.ind4integ = [slice(self.iCompute[d].start,
+                                    stops[d]) for d in xrange(self._dim)]
+            start = [self.position_in_parent[d].start - self.substart[d]
+                     for d in xrange(self._dim)]
+            ## position of the LOCAL submesh in the global grid
+            ## of the submesh (not the grid of the parent mesh!)
+            self.position = [slice(start[d], start[d] + self.resolution[d])
+                             for d in xrange(self._dim)]
+        else:
+            self.position_in_parent = None
+            self.position = None
+            self.iCompute = [slice(0, 0), ] * self._dim
+            self.ind4integ = self.iCompute
+            self.resolution = [0, ] * self._dim
+
+        ## Coordinates of the points of the local submesh
+        self.coords = [None, ] * self._dim
+        self.coords[0] = self.mesh.coords[0][self.iCompute[0]]
+        if self._dim > 1:
+            self.coords[1] = self.mesh.coords[1][:, self.iCompute[1], ...]
+        if self._dim > 2:
+            self.coords[2] = self.mesh.coords[2][:, :, self.iCompute[2]]
+        self.coords = tuple(self.coords)
+        ## Same as self.coords but reduced for integration on the submesh
+        self.coords4int = [None, ] * self._dim
+        self.coords4int[0] = self.mesh.coords[0][self.ind4integ[0]]
+        if self._dim > 1:
+            self.coords4int[1] = self.mesh.coords[1][:, self.ind4integ[1], ...]
+        if self._dim > 2:
+            self.coords4int[2] = self.mesh.coords[2][:, :, self.ind4integ[2]]
+        self.coords4int = tuple(self.coords4int)
+
+    def check_boundaries(self):
+        """
+        Special care when some boundaries of the submesh are on the
+        upper boundaries of the parent mesh.
+        Remind that for periodic bc, such a point does not really
+        exists in the parent mesh.
+        """
+        # List of directions which are periodic
+        periodic_dir = self.mesh.domain.i_periodic_boundaries()
+        if len(periodic_dir) > 0:
+            ll = np.where(
+                self.subend ==
+                self.mesh.global_indices(self.mesh.domain.end))[0]
+            return (ll != self._n_dir).all()
+
+        return True
+
+    def global_resolution(self):
+        """
+        @return the resolution of the global grid (on the whole
+        domain, whatever the mpi distribution is).
+        """
+        return self.discretization.resolution
+
+    def cx(self):
+        return self.coords[0][:, ...]
+
+    def cy(self):
+        assert self._dim > 1
+        return self.coords[1][0, :, ...]
+
+    def cz(self):
+        assert self._dim > 2
+        return self.coords[2][0, 0, :]

HySoP/hysop/domain/tests/test_control_box.py

+from parmepy.domain.subsets.control_box import ControlBox
+from parmepy.domain.subsets.boxes import SubBox
+from parmepy.tools.parameters import Discretization
+from parmepy import Field, Box
+import numpy as np
+import math
+import parmepy.tools.numpywrappers as npw
+
+
+Nx = 128
+Ny = 96
+Nz = 102
+g = 2
+
+
+def v3d(res, x, y, z, t):
+    res[0][...] = 1.
+    res[1][...] = 1.
+    res[2][...] = 1.
+    return res
+
+
+def v2d(res, x, y, t):
+    res[0][...] = 1.
+    res[1][...] = 1.
+    return res
+
+
+def f_test(x, y, z):
+    return 1
+
+
+def f_test_2(x, y, z):
+    return math.cos(z)
+
+ldef = [1.1, 0.76, 1.0]
+discr3D = Discretization([Nx + 1, Ny + 1, Nz + 1], [g - 1, g - 2, g])
+discr2D = Discretization([Nx + 1, Ny + 1], [g - 1, g])
+orig = np.asarray([0.1, -0.3, 0.5])
+ldom = [math.pi * 2., math.pi * 2., math.pi * 2.]
+xdef = orig + 0.2
+
+
+def check_control_box(discr, xr, lr):
+    dim = len(discr.resolution)
+    dom = Box(dimension=dim, length=ldom[:dim],
+              origin=orig[:dim])
+    # Starting point and length of the subdomain
+    cb = ControlBox(origin=xr, length=lr, parent=dom)
+    assert (cb.length == lr).all()
+    assert (cb.origin == xr).all()
+    # discretization of the dom and of its subset
+    topo = dom.create_topology(discr)
+    cb.discretize(topo)
+    assert cb.mesh.values()[0].mesh == topo.mesh
+    assert cb.mesh.keys()[0] == topo
+    assert len(cb.surf) == dim * 2
+    for s in cb.surf:
+        assert isinstance(s, SubBox)
+        ll = s.real_length[topo]
+        assert len(np.where(ll == 0.0)) == 1
+    return topo, cb
+
+
+def test_cb_3D():
+    topo, cb = check_control_box(discr3D, xdef, ldef)
+    velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    i0 = cb.integrate_field_allc(velo, topo)
+    vref = np.prod(cb.real_length[topo])
+    assert (np.abs(i0 - vref) < 1e-6).all()
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref = np.prod(cb.real_length[topo][ilist])
+        isurf = cb.integrate_on_faces(velo, topo, [2 * i])
+        assert np.abs(isurf - sref) < 1e-6
+        isurf = cb.integrate_on_faces(velo, topo, [2 * i + 1])
+        assert np.abs(isurf - sref) < 1e-6
+
+
+def test_cb_3D_2():
+    topo, cb = check_control_box(discr3D, xdef, ldef)
+    velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    list_dir = np.arange(2 * nbc)
+    sref = 0
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref += np.prod(cb.real_length[topo][ilist])
+    sref *= 2.
+    isurf = cb.integrate_on_faces(velo, topo, list_dir)
+    assert np.abs(isurf - sref) < 1e-6
+
+
+def test_cb_3D_3():
+    topo, cb = check_control_box(discr3D, xdef, ldef)
+    velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    list_dir = np.arange(2 * nbc)
+    sref = 0
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref += np.prod(cb.real_length[topo][ilist])
+    sref *= 2.
+    isurf = cb.integrate_on_faces_allc(velo, topo, list_dir)
+    assert (np.abs(isurf - sref) < 1e-6).all()
+
+
+def test_cb_2D():
+    topo, cb = check_control_box(discr2D, xdef[:2], ldef[:2])
+    velo = Field(domain=topo.domain, isVector=True, formula=v2d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    i0 = cb.integrate_field_allc(velo, topo)
+    vref = np.prod(cb.real_length[topo])
+    assert (np.abs(i0 - vref) < 1e-6).all()
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref = np.prod(cb.real_length[topo][ilist])
+        isurf = cb.integrate_on_faces(velo, topo, [2 * i])
+        assert np.abs(isurf - sref) < 1e-6
+        isurf = cb.integrate_on_faces(velo, topo, [2 * i + 1])
+        assert np.abs(isurf - sref) < 1e-6
+
+
+def test_cb_2D_2():
+    topo, cb = check_control_box(discr2D, xdef[:2], ldef[:2])
+    velo = Field(domain=topo.domain, isVector=True, formula=v2d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    list_dir = np.arange(2 * nbc)
+    sref = 0
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref += np.prod(cb.real_length[topo][ilist])
+    sref *= 2.
+    isurf = cb.integrate_on_faces(velo, topo, list_dir)
+    assert np.abs(isurf - sref) < 1e-6
+
+
+def test_cb_2D_3():
+    topo, cb = check_control_box(discr2D, xdef[:2], ldef[:2])
+    velo = Field(domain=topo.domain, isVector=True, formula=v2d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    list_dir = np.arange(2 * nbc)
+    sref = 0
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref += np.prod(cb.real_length[topo][ilist])
+    sref *= 2.
+    isurf = cb.integrate_on_faces_allc(velo, topo, list_dir)
+    assert (np.abs(isurf - sref) < 1e-6).all()
+
+
+def test_cb_3D_fulldomain():
+    """
+    A cb defined on the whole domain.
+    Integrals on upper surfaces must fail,
+    because of periodic bc.
+    """
+    topo, cb = check_control_box(discr3D, orig, ldom)
+    velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
+    velo.discretize(topo)
+    velo.initialize(topo=topo)
+    i0 = cb.integrate_field_allc(velo, topo)
+    vref = np.prod(cb.real_length[topo])
+    assert (np.abs(i0 - vref) < 1e-6).all()
+    nbc = velo.nbComponents
+    sref = npw.zeros(nbc)
+    dirs = np.arange(nbc)
+    for i in xrange(nbc):
+        ilist = np.where(dirs != i)[0]
+        sref = np.prod(cb.real_length[topo][ilist])
+        isurf = cb.integrate_on_faces(velo, topo, [2 * i])
+        assert np.abs(isurf - sref) < 1e-6
+        res = False
+        try:
+            isurf = cb.integrate_on_faces(velo, topo, [2 * i + 1])
+        except:
+            res = True
+        assert res
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_cb_3D()
+    test_cb_3D_2()
+    test_cb_3D_3()
+    test_cb_2D()
+    test_cb_2D_2()
+    test_cb_2D_3()
+    test_cb_3D_fulldomain()