From 10b8ece7918702b94ffdfe1239b22cb9508eb28d Mon Sep 17 00:00:00 2001
From: Keck Jean-Baptiste <jean-baptiste.keck@imag.fr>
Date: Tue, 12 Sep 2017 23:59:08 +0200
Subject: [PATCH] fixed exchange ghosts
---
.gitignore | 1 +
.../taylor_green/taylor_green_monoscale.py | 17 ++--
.../device/codegen/structs/mesh_info.py | 4 +-
.../operator/directional/advection_dir.py | 29 ++++--
hysop/core/mpi/redistribute.py | 4 +-
hysop/fields/cartesian_discrete_field.py | 97 ++++++++++++++++---
hysop/mesh/cartesian_mesh.py | 78 ++++++++++-----
hysop/operator/adapt_timestep.py | 19 ++--
hysop/operator/hdf_io.py | 7 +-
hysop/problem.py | 4 +-
hysop/tools/io_utils.py | 4 +-
hysop/topology/cartesian_topology.py | 7 ++
12 files changed, 202 insertions(+), 69 deletions(-)
diff --git a/.gitignore b/.gitignore
index 10c2abc5b..b86474991 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
*.pyc
+*.pyo
*.so
*.log
hysop.egg-info/
diff --git a/examples/taylor_green/taylor_green_monoscale.py b/examples/taylor_green/taylor_green_monoscale.py
index 820ffed57..1424e7bd1 100644
--- a/examples/taylor_green/taylor_green_monoscale.py
+++ b/examples/taylor_green/taylor_green_monoscale.py
@@ -29,7 +29,7 @@ def init_velocity(data, coords, dim):
elif dim==2:
(x,y) = coords
data[0][...] = +1.0
- data[1][...] = -1.0
+ data[1][...] = -0.0
else:
raise NotImplementedError('dimension {}'.format(dim))
@@ -43,11 +43,11 @@ def init_vorticity(data, coords, dim):
data[2][...] = 2. * sin(x) * sin(y) * cos(z)
elif dim==2:
(x,y) = coords
- data[0][...] = cos(x)*sin(y)
+ data[0][...] = cos(x)*cos(y)
else:
raise NotImplementedError('dimension {}'.format(dim))
-def do_solve(npts, viscosity=1./1600., lcfl=0.125, cfl=2.0):
+def do_solve(npts, viscosity=1./1600., lcfl=0.125, cfl=0.5):
autotuner_config = AutotunerConfig(autotuner_flag=AutotunerFlags.ESTIMATE,
prune_threshold=1.2, override_cache=True, verbose=0)
kernel_config = KernelConfig(autotuner_config=autotuner_config)
@@ -68,7 +68,6 @@ def do_solve(npts, viscosity=1./1600., lcfl=0.125, cfl=2.0):
variables = {velo: d3d, vorti: d3d},
method = {Remesh: Remesh.L2_1},
)
- advec.dump_outputs(directions=1, name_postfix='_X', kill=True)
if dim==3:
stretch = DirectionalStretching(
name='stretch',
@@ -91,22 +90,20 @@ def do_solve(npts, viscosity=1./1600., lcfl=0.125, cfl=2.0):
method={KernelConfig: kernel_config})
splitting.push_operators(advec)#, stretch)
- simu = Simulation(start=0.0, end=10.0, dt0=2*pi/8.0)
+ simu = Simulation(start=0.0, end=10.0, dt0=0.01)
adapt_dt = AdaptiveTimeStep(simu.dt, dt_coeff=0.95)
adapt_dt.push_cfl_criteria(cfl, velo)
adapt_dt.push_lcfl_criteria(lcfl, vorti)
problem = Problem(method={ComputeGranularity:0})
- problem.insert(splitting, diffusion, poisson, adapt_dt)
+ problem.insert(splitting, adapt_dt)
problem.dump_inputs()
problem.build()
- problem.display()
+ #problem.display()
dfields = problem.get_input_discrete_fields()
dfields[velo].initialize(formula=init_velocity, dim=dim)
dfields[vorti].initialize(formula=init_vorticity, dim=dim)
-
- print 'DX=', dfields[velo].mesh.space_step
with printoptions(threshold=10000, linewidth=240,
nanstr='nan', infstr='inf',
@@ -115,6 +112,6 @@ def do_solve(npts, viscosity=1./1600., lcfl=0.125, cfl=2.0):
problem.finalize()
if __name__=='__main__':
- N = 8 + 1
+ N = 64 + 1
do_solve(N)
diff --git a/hysop/backend/device/codegen/structs/mesh_info.py b/hysop/backend/device/codegen/structs/mesh_info.py
index d7473278f..a04104a3d 100644
--- a/hysop/backend/device/codegen/structs/mesh_info.py
+++ b/hysop/backend/device/codegen/structs/mesh_info.py
@@ -229,8 +229,8 @@ class MeshInfoStruct(StructCodeGenerator):
dx = mesh.space_step[::-1]
ghosts = mesh.ghosts[::-1]
- gresolution = mesh.global_resolution[::-1]
- gcompute_resolution = mesh.global_resolution[::-1]
+ gresolution = mesh.grid_resolution[::-1]
+ gcompute_resolution = mesh.grid_resolution[::-1]
gxmin = mesh.global_origin[::-1]
gxmax = mesh.global_end[::-1]
gsize = mesh.global_length[::-1]
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index 05bb27aaf..240148e3b 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -90,7 +90,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
lidx[...] %= N
elif __debug__:
min_lidx, max_lidx = lidx.min(), lidx.max()
- M = Vin.shape[-1]-1
+ M = Vin.shape[-1]
assert min_lidx >= 0, 'index out of bounds: {} < 0'.format(min_lidx)
assert max_lidx < M-1, 'index out of bounds: {} >= {}'.format(max_lidx, M-1)
ridx[...] = lidx
@@ -222,9 +222,9 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
N = R0.shape[-1]
# print 'IDX=\n{}'.format('\n'.join([str(x[0]) for x in tuple(mesh.iter_compute_axes(cr))]))
- print 'INPUT'
- print sinputs.values()[0][0][sinputs.values()[0].compute_slices]
- print
+ #print 'INPUT'
+ #print sinputs.values()[0][0][sinputs.values()[0].compute_slices]
+ #print
for (idx,_,_,_) in mesh.iter_compute_axes(cr):
Pi = pos[idx].lock()
R0[...] = Pi
@@ -291,8 +291,25 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
sout[...] = Si
# print 'SOUT \n', sout
sid+=1
- print 'FINAL'
- print soutputs.values()[0][0][soutputs.values()[0].compute_slices]
+
+ # exchange and sum remeshed ghosts
+ for sout in soutputs.values():
+ buffers = sout.buffers
+ outer_ghosts = sout.outer_ghost_slices
+ inner_ghosts = sout.inner_ghost_slices
+ nprocs = sout.mesh.proc_shape
+
+ for buf in buffers:
+ lsrc, rsrc = inner_ghosts[-1]
+ ldst, rdst = outer_ghosts[-1]
+ if (nprocs[-1] == 1):
+ buf[lsrc] += buf[rdst]
+ buf[rsrc] += buf[ldst]
+ buf[ldst] = buf[rsrc]
+ buf[rdst] = buf[lsrc]
+ else:
+ msg='MPI exchange not implemented yet.'
+ raise RuntimeError(msg)
@debug
diff --git a/hysop/core/mpi/redistribute.py b/hysop/core/mpi/redistribute.py
index b940c84bb..cfdaa2dea 100644
--- a/hysop/core/mpi/redistribute.py
+++ b/hysop/core/mpi/redistribute.py
@@ -150,8 +150,8 @@ class RedistributeIntra(RedistributeOperatorBase):
return False
# source and target must have the same global resolution
- source_res = tin.mesh.global_resolution
- target_res = tout.mesh.global_resolution
+ source_res = tin.mesh.grid_resolution
+ target_res = tout.mesh.grid_resolution
if not npw.array_equal(source_res, target_res):
return False
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index 999ca562c..1df5c591f 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -84,10 +84,6 @@ class CartesianDiscreteFieldView(DiscreteFieldView):
for i,var in enumerate(data):
if backend.can_wrap(var):
data[i] = backend.wrap(var)
- # else:
- # msg='Unknown data type {} for backend {}.'
- # msg=msg.format(var.__class__, backend.kind)
- # raise ValueError(msg)
for idx, _data in zip(components, data):
if _data.ndim != self.dim:
@@ -212,9 +208,19 @@ class CartesianDiscreteFieldView(DiscreteFieldView):
def _get_compute_slices(self):
"""Return a tuple of slices indexing the local compute mesh."""
return self.mesh.local_compute_slices
- def _get_ghost_slices(self, directions=None, ):
- """Return a tuple of tuples of slices indexing the ghosts."""
- return self.mesh.local_ghost_slices
+ def _get_outer_ghost_slices(self):
+ """
+ Return a tuple of tuples of slices indexing the local outer ghosts.
+ Those slices corresponds to local to process ghosts.
+ """
+ return self.mesh.local_outer_ghost_slices
+ def _get_inner_ghost_slices(self):
+ """
+ Return a tuple of tuples of slices indexing the local inner ghosts.
+ Those slices corresponds to neighbour processes overlap on local
+ compute slices.
+ """
+ return self.mesh.local_inner_ghost_slices
def short_description(self):
"""Short description of this discrete field."""
@@ -308,6 +314,9 @@ CartesianDiscreteFieldView (id={}, tag={})
# call formula
formula(**formula_kwds)
+ # update ghosts
+ self.exchange_ghosts(data=data)
+
assert all((d.size == self.size) for d in data), 'Array size was altered.'
assert all((d.dtype == self.dtype) for d in data), 'Array dtype was altered.'
assert all(all(d.shape == self.shape) for d in data), 'Array shape was altered.'
@@ -334,6 +343,70 @@ CartesianDiscreteFieldView (id={}, tag={})
self.topology.cart_comm.Allreduce(result, gresult)
return gresult ** (1.0/p)
+ def print_with_ghosts(self, component=0, compute=None,
+ inner_ghosts=None, outer_ghosts='X', **print_opts):
+ """
+ Print values with ghosts replaced by symbol.
+ Mainly for debug purposes.
+ """
+ data = self.data[component].get()
+ strarr = np.empty_like(data, dtype=object)
+ strarr[...] = data
+ if (compute is not None):
+ strarr[self.compute_slices] = compute
+ if (inner_ghosts is not None):
+ for lg, rg in self.inner_ghost_slices:
+ strarr[lg] = inner_ghosts
+ strarr[rg] = inner_ghosts
+ if (outer_ghosts is not None):
+ for lg, rg in self.outer_ghost_slices:
+ strarr[lg] = outer_ghosts
+ strarr[rg] = outer_ghosts
+
+ _formatter = {
+ object: lambda x: '{:^6}'.format(x)[:6],
+ float: lambda x: '{:6.2f}'.format(x)
+ }
+
+ _print_opts = dict(threshold=10000, linewidth=1000,
+ nanstr='nan', infstr='inf',
+ formatter={'object': lambda x: _formatter.get(type(x),_formatter[object])(x)})
+ _print_opts.update(print_opts)
+
+ from hysop.tools.contexts import printoptions
+ with printoptions(**_print_opts):
+ print strarr
+
+
+ def exchange_ghosts(self, components=None, directions=None, data=None,
+ custom_op=None):
+ """Exchange ghosts on specified components and directions."""
+ from hysop.core.arrays.array import Array
+ components = to_tuple(components or range(self.nb_components))
+ directions = to_tuple(directions or range(self.dim))
+
+ data = data or self.data
+ data = tuple(d.handle if isinstance(d,Array) else d for d in data)
+ # exchange_ghosts currently only supports hosts buffers
+ check_instance(data, tuple, values=np.ndarray)
+
+ outer_ghosts = self.outer_ghost_slices
+ inner_ghosts = self.inner_ghost_slices
+ nprocs = self.mesh.proc_shape
+
+ for component in components:
+ field = self.data[component]
+ for direction in directions:
+ lsrc, rsrc = inner_ghosts[direction]
+ ldst, rdst = outer_ghosts[direction]
+ if (nprocs[direction] == 1):
+ field[rdst] = field[lsrc]
+ field[ldst] = field[rsrc]
+ else:
+ msg='MPI exchange not implemented yet.'
+ raise RuntimeError(msg)
+
+
size = property(_get_size)
shape = property(_get_shape)
compute_resolution = property(_get_compute_resolution)
@@ -341,7 +414,8 @@ CartesianDiscreteFieldView (id={}, tag={})
ghosts = property(_get_ghosts)
compute_slices = property(_get_compute_slices)
- ghost_slices = property(_get_ghost_slices)
+ inner_ghost_slices = property(_get_inner_ghost_slices)
+ outer_ghost_slices = property(_get_outer_ghost_slices)
axes = property(_get_axes)
tstate = property(_get_tstate)
@@ -352,6 +426,7 @@ CartesianDiscreteFieldView (id={}, tag={})
+
class CartesianDiscreteField(CartesianDiscreteFieldView, DiscreteField):
"""
Discrete representation of cartesian scalar or vector fields,
@@ -431,10 +506,8 @@ class CartesianDiscreteField(CartesianDiscreteFieldView, DiscreteField):
if (vd is not None):
data[i][topology.mesh.local_compute_slices] = vd
if (vg is not None):
- for ls,rs in topology.mesh.local_ghost_slices:
- print ls, vg
- data[i][ls] = 0
- # data[i][rs] = vg
+ for ls,rs in topology.mesh.outer_ghost_slices:
+ data[i][rs] = vg
if self.topology_state.is_read_only:
npw.set_readonly(data[i])
diff --git a/hysop/mesh/cartesian_mesh.py b/hysop/mesh/cartesian_mesh.py
index 1f50fdd88..a8479502f 100644
--- a/hysop/mesh/cartesian_mesh.py
+++ b/hysop/mesh/cartesian_mesh.py
@@ -64,9 +64,17 @@ class CartesianMeshView(MeshView):
"""Return True if the current process has points on the current mesh."""
return self._mesh._on_proc
+ def _get_grid_resolution(self):
+ """
+ Effective resolution of the global mesh.
+ Corresponds to self.global_resolution - topology.is_periodic.
+ See global_resolution for topology global resolution.
+ """
+ return self.__get_transposed_mesh_attr('_grid_resolution')
def _get_global_resolution(self):
"""
- Global resolution of the mesh.
+ Global resolution of the associated topology.
+ See grid_resolution for effective global mesh size.
"""
return self.__get_transposed_mesh_attr('_global_resolution')
def _get_global_start(self):
@@ -197,14 +205,23 @@ class CartesianMeshView(MeshView):
This is the local computational view (ie. local indices of the computational points).
"""
return self.__get_transposed_mesh_attr('_local_compute_slices')
- def _get_local_ghost_slices(self):
+ def _get_local_inner_ghost_slices(self):
"""
Return a list of slices defining the ghosts in this arrays as local indices.
- There are two ghost slices per specified directions.
+ Those slices corresponds to neighbour processes overlap on local
+ compute slices.
For each direction, a pair of tuples of slices is returned, one
- for left and right ghosts.
+ for left and one for the right ghosts.
"""
- return self.__get_transposed_mesh_attr('_local_ghost_slices')
+ return self.__get_transposed_mesh_attr('_local_inner_ghost_slices')
+ def _get_local_outer_ghost_slices(self):
+ """
+ Return a list of slices defining the ghosts in this arrays as local indices.
+ Those slices corresponds to local to process ghosts.
+ For each direction, a pair of tuples of slices is returned, one
+ for left and one for the right ghosts.
+ """
+ return self.__get_transposed_mesh_attr('_local_outer_ghost_slices')
def _get_local_origin(self):
"""
Origin of the first computational point on the physical box
@@ -452,6 +469,8 @@ class CartesianMeshView(MeshView):
topology = property(_get_topology)
on_proc = property(_get_on_proc)
+ grid_resolution=property(_get_grid_resolution)
+
global_resolution=property(_get_global_resolution)
global_start=property(_get_global_start)
global_stop=property(_get_global_stop)
@@ -466,7 +485,8 @@ class CartesianMeshView(MeshView):
local_start=property(_get_local_start)
local_stop=property(_get_local_stop)
local_compute_slices=property(_get_local_compute_slices)
- local_ghost_slices=property(_get_local_ghost_slices)
+ local_inner_ghost_slices=property(_get_local_inner_ghost_slices)
+ local_outer_ghost_slices=property(_get_local_outer_ghost_slices)
local_origin=property(_get_local_origin)
local_end=property(_get_local_end)
local_length=property(_get_local_length)
@@ -602,6 +622,9 @@ class CartesianMesh(CartesianMeshView, Mesh):
global_start = npw.asintegerarray(global_start).copy()
assert global_start.size == dim
assert (global_start>=0).all()
+
+ # Remove 1 point on each periodic axe because of periodicity
+ grid_resolution = global_resolution - domain.periodicity
# ---
# Attributes relative to the distributed mesh
@@ -610,7 +633,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
local_resolution = npw.asintegerarray(local_resolution)
assert local_resolution.size == dim
assert (local_resolution > 2*ghosts).all()
- assert (local_resolution <= global_resolution + 2*ghosts).all()
+ assert (local_resolution <= grid_resolution + 2*ghosts).all()
compute_resolution = local_resolution - 2*ghosts
# deduce all other attributes
@@ -636,14 +659,21 @@ class CartesianMesh(CartesianMeshView, Mesh):
local_start = ghosts
local_stop = local_resolution - ghosts
local_compute_slices = tuple( slice(i,j) for (i,j) in zip(local_start, local_stop) )
- local_ghost_slices = []
+ local_inner_ghost_slices = []
+ local_outer_ghost_slices = []
for i in xrange(dim):
- lslices = tuple( slice(None) if j!=i else slice(local_start[i]-ghosts[i],
+ inner_lslices = tuple( slice(None) if j!=i else slice(local_start[i],
+ local_start[i]+ghosts[i]) for j in xrange(dim) )
+ inner_rslices = tuple( slice(None) if j!=i else slice(local_stop[i]-ghosts[i],
+ local_stop[i]) for j in xrange(dim) )
+ outer_lslices = tuple( slice(None) if j!=i else slice(local_start[i]-ghosts[i],
local_start[i]) for j in xrange(dim) )
- rslices = tuple( slice(None) if j!=i else slice(local_stop[i],
+ outer_rslices = tuple( slice(None) if j!=i else slice(local_stop[i],
local_stop[i]+ghosts[i]) for j in xrange(dim) )
- local_ghost_slices.append( (lslices,rslices) )
- local_ghost_slices = tuple(local_ghost_slices)
+ local_inner_ghost_slices.append( (inner_lslices, inner_rslices) )
+ local_outer_ghost_slices.append( (outer_lslices, outer_rslices) )
+ local_inner_ghost_slices = tuple(local_inner_ghost_slices)
+ local_outer_ghost_slices = tuple(local_outer_ghost_slices)
local_indices = tuple(np.arange(Ni, dtype=HYSOP_INTEGER) for Ni in local_resolution)
local_coords = tuple(
@@ -670,15 +700,13 @@ class CartesianMesh(CartesianMeshView, Mesh):
rboundaries = np.asarray( [bc if at_right else BoundaryCondition.NONE
for (bc, at_right) in zip(global_boundaries[1], is_at_right_boundary) ])
local_boundaries = (lboundaries, rboundaries)
- is_periodic = np.logical_and(lboundaries==rboundaries,
- lboundaries==BoundaryCondition.PERIODIC)
- npw.set_readonly(global_resolution, global_start, global_stop,
+ npw.set_readonly(global_resolution, grid_resolution, global_start, global_stop,
global_origin, global_end, global_length,
global_boundaries[0], global_boundaries[1],
local_resolution, local_start, local_stop,
local_origin, local_end, local_length,
- local_boundaries[0], local_boundaries[1], is_periodic,
+ local_boundaries[0], local_boundaries[1],
compute_resolution, ghosts,
proc_coords, proc_shape,
is_at_left_boundary, is_at_right_boundary,
@@ -688,6 +716,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
npw.set_readonly(*local_compute_indices)
npw.set_readonly(*local_compute_coords)
+ self._grid_resolution = grid_resolution
self._global_resolution = global_resolution
self._global_start = global_start
self._global_stop = global_stop
@@ -701,8 +730,6 @@ class CartesianMesh(CartesianMeshView, Mesh):
self._local_resolution = local_resolution
self._local_start = local_start
self._local_stop = local_stop
- self._local_compute_slices = local_compute_slices
- self._local_ghost_slices = local_ghost_slices
self._local_origin = local_origin
self._local_end = local_end
self._local_length = local_length
@@ -712,6 +739,10 @@ class CartesianMesh(CartesianMeshView, Mesh):
self._local_compute_indices = local_compute_indices
self._local_compute_coords = local_compute_coords
+ self._local_compute_slices = local_compute_slices
+ self._local_inner_ghost_slices = local_inner_ghost_slices
+ self._local_outer_ghost_slices = local_outer_ghost_slices
+
self._compute_resolution = compute_resolution
self._ghosts = ghosts
self._proc_coords = proc_coords
@@ -719,7 +750,6 @@ class CartesianMesh(CartesianMeshView, Mesh):
self._space_step = space_step
self._is_at_left_boundary = is_at_left_boundary
self._is_at_right_boundary = is_at_right_boundary
- self._is_periodic = is_periodic
self._on_proc = True
if __debug__:
@@ -729,6 +759,7 @@ class CartesianMesh(CartesianMeshView, Mesh):
# check variables and properties at the same time
from hysop.topology.cartesian_topology import CartesianTopologyView
check_instance(self.topology, CartesianTopologyView)
+ check_instance(self.grid_resolution, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
check_instance(self.global_resolution, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
check_instance(self.global_start, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
check_instance(self.global_stop, np.ndarray, dtype=HYSOP_INTEGER, shape=(dim,))
@@ -751,8 +782,12 @@ class CartesianMesh(CartesianMeshView, Mesh):
check_instance(self.local_end, np.ndarray, dtype=HYSOP_REAL, shape=(dim,))
check_instance(self.local_length, np.ndarray, dtype=HYSOP_REAL, shape=(dim,))
check_instance(self.local_compute_slices, tuple, values=slice, size=dim)
- check_instance(self.local_ghost_slices, tuple, values=tuple, size=dim)
- for lg,rg in self.local_ghost_slices:
+ check_instance(self.local_inner_ghost_slices, tuple, values=tuple, size=dim)
+ for lg,rg in self.local_inner_ghost_slices:
+ check_instance(lg, tuple, values=slice, size=dim)
+ check_instance(rg, tuple, values=slice, size=dim)
+ check_instance(self.local_outer_ghost_slices, tuple, values=tuple, size=dim)
+ for lg,rg in self.local_outer_ghost_slices:
check_instance(lg, tuple, values=slice, size=dim)
check_instance(rg, tuple, values=slice, size=dim)
check_instance(self.local_boundaries, tuple, values=np.ndarray, size=2)
@@ -794,7 +829,6 @@ class CartesianMesh(CartesianMeshView, Mesh):
check_instance(self.is_at_left_boundary, np.ndarray, dtype=bool, shape=(dim,))
check_instance(self.is_at_right_boundary, np.ndarray, dtype=bool, shape=(dim,))
check_instance(self.is_at_boundary, np.ndarray, dtype=bool, shape=(dim,))
- # check_instance(self.is_periodic, np.ndarray, dtype=bool, shape=(dim,))
check_instance(self.on_proc, bool)
def view(self, topology_state):
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index e41959748..fa218ecd5 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -107,22 +107,25 @@ class CflTimestepCriteria(TimestepCriteria):
def compute_criteria(self, **kwds):
from hysop.fields.cartesian_discrete_field import CartesianDiscreteFieldView
+
dv = self.dvelocity
cfl = self.cfl
+ dt = 0.0
+
if isinstance(dv, CartesianDiscreteFieldView):
dx = dv.dfield.mesh.space_step
-
- V0_inf = npw.maximum( npw.abs(dv[0].max()), npw.abs(dv[0].min()) )
- dt = (V0_inf*dx[0]) / cfl
- for i in xrange(dv.nb_components-1):
- Vi_inf = npw.maximum( npw.abs(dv[i].max()), npw.abs(dv[i].min()) )
- dti = (Vi_inf*dx[i]) / cfl
- dt = npw.maximum(dt, dti)
- return dt
+ for i in xrange(dv.nb_components):
+ Vi_min, Vi_max = dv[i].min(), dv[i].max()
+ Vi_inf = max(abs(Vi_min), abs(Vi_max))
+ dti = (dx[i]*cfl) / Vi_inf
+ dt = max(dt, dti)
else:
msg='CFL criteria has not been implemented for {} discrete fields.'
msg=msg.format(type(self.dvelocity))
raise RuntimeError(msg)
+
+ assert dt>0.0
+ return dt
class MergeTimeStepCriterias(TimestepCriteria):
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index 938654319..f40541064 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -8,7 +8,7 @@
"""
from abc import ABCMeta, abstractmethod
-from hysop.deps import h5py
+from hysop.deps import h5py, sys
from hysop.core.graph.graph import discretized
from hysop.constants import DirectionLabels, HYSOP_REAL
from hysop.tools.decorators import debug
@@ -130,7 +130,7 @@ class HDF_IO(ComputationalGraphOperator):
# Global resolution for hdf5 output (warning : this must
# be the whole domain resolution, not the subset resolution)
- self._global_resolution = list(refmesh.global_resolution)
+ self._global_resolution = list(refmesh.grid_resolution)
self._slices = refmesh.local_compute_slices
if refmesh.on_proc:
sl = list(refmesh.global_compute_slices)
@@ -235,11 +235,10 @@ class HDF_Writer(HDF_IO):
if simulation is None:
raise ValueError("Missing simulation value for monitoring.")
ite = simulation.current_iteration
- if ite == -1 or ite % self.io_params.frequency == 0:
+ if (ite == -1) or (ite % self.io_params.frequency == 0):
self.step(simulation)
self._count += 1
if self.kill:
- import sys
sys.exit(1)
def createXMFFile(self):
diff --git a/hysop/problem.py b/hysop/problem.py
index 45b9df2fc..4cd748dc0 100644
--- a/hysop/problem.py
+++ b/hysop/problem.py
@@ -26,7 +26,7 @@ class Problem(ComputationalGraph):
self.setup(work)
@debug
- def solve(self, simu, **kargs):
+ def solve(self, simu, max_iter=None, **kargs):
vprint('\nSolving problem...')
simu.initialize()
while not simu.is_over:
@@ -34,6 +34,8 @@ class Problem(ComputationalGraph):
simu.print_state()
self.apply(simulation=simu, **kargs)
simu.advance()
+ if (max_iter is not None) and (simu.current_iteration > max_iter):
+ break
simu.finalize()
@debug
diff --git a/hysop/tools/io_utils.py b/hysop/tools/io_utils.py
index 83292409e..dcdc35af4 100755
--- a/hysop/tools/io_utils.py
+++ b/hysop/tools/io_utils.py
@@ -461,10 +461,10 @@ class XMF(object):
# Check substet to find the required grid resolution
dx = list(topo.mesh.space_step)
if subset is not None:
- res = list(subset.mesh[topo].global_resolution)
+ res = list(subset.mesh[topo].grid_resolution)
orig = list(subset.real_orig[topo])
else:
- res = list(topo.mesh.global_resolution)
+ res = list(topo.mesh.grid_resolution)
orig = list(topo.domain.origin)
resolution = [1,]*3
origin = [0.0,]*3
diff --git a/hysop/topology/cartesian_topology.py b/hysop/topology/cartesian_topology.py
index f022fcba8..4ead103ad 100644
--- a/hysop/topology/cartesian_topology.py
+++ b/hysop/topology/cartesian_topology.py
@@ -880,6 +880,13 @@ class CartesianTopology(CartesianTopologyView, Topology):
local_resolution = nbpoints.copy()
local_resolution += 2 * discretization.ghosts
+
+ msg='\nLocal compute shape is smaller than the total number of ghosts, '
+ msg+='on at least one axis:'
+ msg+=' *compute shape: {}'.format(nbpoints)
+ msg+=' *ghosts: 2*{}'.format(discretization.ghosts)
+ if not np.all(nbpoints >= 2 * discretization.ghosts):
+ raise RuntimeError(msg)
# Global indices for the local mesh points
global_start = proc_coords * pts_noghost
--
GitLab