diff --git a/HySoP/hysop/domain/mesh.py b/HySoP/hysop/domain/mesh.py
index 419cd4ed7008b69d3b5cfde855d70494358a6393..2332da396e4340add5694e7186448ca65cbae2a3 100644
--- a/HySoP/hysop/domain/mesh.py
+++ b/HySoP/hysop/domain/mesh.py
@@ -273,11 +273,18 @@ class SubMesh(object):
## position of the submesh in the global grid of its parent mesh.
self.global_position_in_parent = [slice(substart[d], subend[d] + 1)
for d in xrange(self._dim)]
+ ## Coordinates of the lowest point of this submesh
+ hh = self.mesh.space_step
+ 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
gres = np.minimum(gres, self.mesh.discretization.resolution - 1)
+ self._s_dir = np.where(gres == 1)[0]
## discretization of the subset
- self.discretization = Discretization(gres, mesh.discretization.ghosts)
+ ## Warning : no ghosts on the submesh!
+ self.discretization = Discretization(gres)
dimension = len(self.substart)
# Find which part of submesh is on the current process and
# find its computational points. Warning:
@@ -307,9 +314,9 @@ class SubMesh(object):
stop[is_last] -= 1
## Same as self.iCompute but recomputed to be used
## for integration on the submesh
+ stop[self._s_dir] += 1
self.ind4integ = [slice(self.iCompute[d].start,
- max(stop[d], self.iCompute[d].start + 1))
- for d in xrange(dimension)]
+ stop[d]) for d in xrange(dimension)]
start = [self.position_in_parent[d].start - self.substart[d]
for d in xrange(dimension)]
@@ -328,7 +335,7 @@ class SubMesh(object):
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], :]
+ 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)
@@ -336,7 +343,7 @@ class SubMesh(object):
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], :]
+ 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)
@@ -348,4 +355,13 @@ class SubMesh(object):
"""
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, :]
diff --git a/HySoP/hysop/domain/subsets/boxes.py b/HySoP/hysop/domain/subsets/boxes.py
index 06195693f5c1293728e1f47753dd22b4c0048e07..244a86b7097930ed8cc3fd2931e604b6bd08bee0 100644
--- a/HySoP/hysop/domain/subsets/boxes.py
+++ b/HySoP/hysop/domain/subsets/boxes.py
@@ -27,7 +27,6 @@ class SubBox(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
@@ -52,9 +51,8 @@ class SubBox(Subset):
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
+ self.real_length[topo] = self.mesh[topo].global_length
+ self.real_orig[topo] = self.mesh[topo].global_origin
def apply(self, tab):
pass
diff --git a/HySoP/hysop/domain/tests/test_submesh.py b/HySoP/hysop/domain/tests/test_submesh.py
index f3135b6ebfe7a657c066162025a5a241a4106bb1..3b2060f90d268ff28d8d8dd5a769dd17cc17bb5d 100644
--- a/HySoP/hysop/domain/tests/test_submesh.py
+++ b/HySoP/hysop/domain/tests/test_submesh.py
@@ -6,48 +6,91 @@ from parmepy.domain.box import Box
from parmepy.tools.parameters import Discretization
from parmepy.domain.mesh import SubMesh
from parmepy.tools.misc import utils
+import numpy as np
Nx = Ny = Nz = 32
g = 2
-def test_submesh():
+def check_submesh(discr):
"""
Periodic mesh
"""
- discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
- gh = discr.ghosts
- dim = gh.size
+ dim = len(discr.resolution)
dom = Box(dimension=dim, origin=[0.1, ] * dim)
- topo = dom.create_topology(discr, cutdir=[False, False, True])
+ topo = dom.create_topology(discr)
mesh = topo.mesh
# Set a starting point and a length for the submesh
- hh = mesh.space_step
- xr = 1. * hh
+ xr = dom.origin + 3. * mesh.space_step
# Find the indices of these points
gstart = mesh.global_indices(xr)
- gend = mesh.global_indices(xr + 15 * hh)
+ gend = mesh.global_indices(xr + 15 * mesh.space_step)
# 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()
+ assert (subm.substart == gstart).all()
+ assert (subm.subend == gend).all()
+ assert np.allclose(subm.global_origin,
+ dom.origin + gstart * mesh.space_step)
# 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()
+ assert (subm.discretization.ghosts == 0).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)]
+ assert subm.position_in_parent == [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
+ pos = subm.position_in_parent
+ subpos = [slice(pos[d].start - gstart[d],
+ pos[d].stop - gstart[d]) for d in xrange(dim)]
+ assert subm.global_position_in_parent == gpos
+ assert subm.position == subpos
+ check_coords(mesh, subm)
+
+
+def check_coords(m1, m2):
+ dim = m1.domain.dimension
+ x0 = np.zeros(dim)
+ xend = np.zeros(dim)
+ x0[0] = m2.cx()[0]
+ xend[0] = m2.cx()[-1]
+ if dim > 1:
+ x0[1] = m2.cy()[0]
+ xend[1] = m2.cy()[-1]
+ if dim > 2:
+ x0[2] = m2.cz()[0]
+ xend[2] = m2.cz()[-1]
+ req_orig = np.maximum(m1.origin + m1.discretization.ghosts * m1.space_step,
+ m1.domain.origin + m2.substart * m1.space_step)
+ req_end = np.minimum(m1.end - m1.discretization.ghosts * m1.space_step,
+ m1.domain.origin + m2.subend * m1.space_step)
+ assert np.allclose(x0, req_orig)
+ assert np.allclose(xend, req_end)
+
+
+def test_submesh_3D():
+ """
+ 3D Periodic mesh
+ """
+ discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
+ check_submesh(discr)
+
+
+def test_submesh_2D():
+ """
+ 2D Periodic mesh
+ """
+ discr = Discretization([Nx + 1, Nz + 1], [g, g])
+ check_submesh(discr)
if __name__ == '__main__':
- test_submesh()
+ test_submesh_3D()
+ test_submesh_2D()
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
index e0a29167942b98910cf18f49db919a7cce1dee78..c7c342ef2681e9c81455f006cbee94a4ffd42cce 100644
--- a/HySoP/hysop/domain/tests/test_subset.py
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -5,7 +5,9 @@ import numpy as np
import math
-Nx = Ny = Nz = 128
+Nx = 128
+Ny = 96
+Nz = 102
g = 2
@@ -16,6 +18,12 @@ def v3d(res, x, y, z, t):
return res
+def v2d(res, x, y, t):
+ res[0][...] = 1.
+ res[1][...] = 1.
+ return res
+
+
def f_test(x, y, z):
return 1
@@ -23,56 +31,49 @@ def f_test(x, y, z):
def f_test_2(x, y, z):
return math.cos(z)
+xdef = [0.2, 0.4, 0.1]
+ldef = [0.3, 0.4, 1.0]
+discr3D = Discretization([Nx + 1, Ny + 1, Nz + 1], [g - 1, g - 2, g])
+discr2D = Discretization([Nx + 1, Ny + 1], [g - 1, g])
-def test_subbox():
- discr = Discretization([Nx + 1, Ny + 1, Nz + 1], [g, g, g])
- dim = 3
- dom = Box(dimension=dim, origin=[0.1, ] * dim)
+
+def check_subset(discr, xr, lr):
+ dim = len(discr.resolution)
+ dom = Box(dimension=dim, length=[math.pi * 2., ] * 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])
+ topo = dom.create_topology(discr)
cub.discretize(topo)
assert cub.mesh.values()[0].mesh == topo.mesh
assert cub.mesh.keys()[0] == topo
+ return topo, cub
+
+
+def test_subbox():
+ check_subset(discr3D, xdef, ldef)
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')
+ topo, cub = check_subset(discr3D, xdef, ldef)
+ velo = Field(domain=topo.domain, 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 False
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])
+ topo, cub = check_subset(discr3D, xdef, ldef)
# 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')
+ #xr = [0.1, 0.2, 0.15]
+ #lr = [0.51, 0.53, 0.3]
+ velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
velo.discretize(topo)
velo.initialize(topo=topo)
vref = np.prod(cub.real_length[topo])
@@ -82,16 +83,10 @@ def test_integ_2():
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])
+ topo, cub = check_subset(discr3D, xdef, ldef)
# 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)
+ #xr = [0.1, 0.2, 0.15]
+ #lr = [0.51, 0.53, 0.3]
nbc = 1
vref = np.prod(cub.real_length[topo])
for i in xrange(1, nbc + 1):
@@ -100,17 +95,10 @@ def test_integ_3():
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])
+ xr = [math.pi * 0.5, ] * 3
+ lr = [math.pi, ] * 3
+ topo, cub = check_subset(discr3D, xr, lr)
# 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):
@@ -119,21 +107,12 @@ def test_integ_4():
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')
+ """
+ subset == plane in 3D domain
+ Test integration on this subset.
+ """
+ topo, cub = check_subset(discr3D, xdef, [0., 0.2, 0.8])
+ velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
velo.discretize(topo)
velo.initialize(topo=topo)
i0 = cub.integrate_field_allc(velo, topo)
@@ -142,21 +121,12 @@ def test_subbox_2d():
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')
+ """
+ subset == line in 3D domain
+ Test integration on this subset.
+ """
+ topo, cub = check_subset(discr3D, xdef, [0., 0., 0.8])
+ velo = Field(domain=topo.domain, isVector=True, formula=v3d, name='velo')
velo.discretize(topo)
velo.initialize(topo=topo)
i0 = cub.integrate_field_allc(velo, topo)
@@ -164,6 +134,31 @@ def test_line():
assert (np.abs(i0 - vref) < 1e-6).all()
+def test_2d_subbox():
+ """
+ subset in 2D domain
+ """
+ topo, cub = check_subset(discr2D, xdef[:2], ldef[:2])
+ velo = Field(domain=topo.domain, isVector=True, formula=v2d, 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_2d_line():
+ """
+ subset == line in 2D domain
+ Test integration on this subset.
+ """
+ topo, cub = check_subset(discr2D, xdef[:2], lr=[0., 1.0])
+ velo = Field(domain=topo.domain, isVector=True, formula=v2d, name='velo')
+ velo.discretize(topo)
+ velo.initialize(topo=topo)
+ i0 = cub.integrate_field_allc(velo, topo)
+ vref = cub.real_length[topo][1]
+ assert (np.abs(i0 - vref) < 1e-6).all()
# This may be useful to run mpi tests
if __name__ == "__main__":
@@ -173,3 +168,5 @@ if __name__ == "__main__":
test_integ_3()
test_subbox_2d()
test_line()
+ test_2d_subbox()
+ test_2d_line()