From 99dc739ff65f3b47859553c367292bd3b03fa75a Mon Sep 17 00:00:00 2001
From: Keck Jean-Baptiste <jean-baptiste.keck@imag.fr>
Date: Sat, 16 Sep 2017 00:43:14 +0200
Subject: [PATCH] io and advec remesh fully optimized
---
.../taylor_green/taylor_green_monoscale.py | 2 +-
hysop/__init__.py | 2 +-
hysop/backend/host/host_buffer.py | 2 +-
.../operator/directional/advection_dir.py | 124 ++++++++++++------
.../backend/host/python/operator/transpose.py | 8 +-
hysop/core/graph/computational_graph.py | 32 +++--
hysop/core/graph/graph.py | 24 ++--
hysop/core/memory/memory_request.py | 6 +-
hysop/fields/cartesian_discrete_field.py | 29 +++-
hysop/mesh/cartesian_mesh.py | 37 ++++--
hysop/operator/adapt_timestep.py | 41 +++---
hysop/operator/hdf_io.py | 94 ++++++++++---
hysop/parameters/numpy_parameter.py | 4 +-
hysop/parameters/scalar_parameter.py | 3 +-
hysop/problem.py | 2 +-
hysop/tools/decorators.py | 4 +-
hysop/tools/io_utils.py | 47 ++-----
17 files changed, 300 insertions(+), 161 deletions(-)
diff --git a/examples/taylor_green/taylor_green_monoscale.py b/examples/taylor_green/taylor_green_monoscale.py
index 42a77492f..0dbde7634 100644
--- a/examples/taylor_green/taylor_green_monoscale.py
+++ b/examples/taylor_green/taylor_green_monoscale.py
@@ -101,7 +101,7 @@ def run(npts=128+1, viscosity=1./1600., lcfl=0.125, cfl=0.5):
problem = Problem(method={ComputeGranularity:0})
problem.insert(splitting, adapt_dt)
- problem.dump_inputs()
+ problem.dump_inputs(fields=vorti)
problem.build()
# problem.display()
diff --git a/hysop/__init__.py b/hysop/__init__.py
index b29be2bef..89f06b50e 100644
--- a/hysop/__init__.py
+++ b/hysop/__init__.py
@@ -20,7 +20,7 @@ __VERBOSE__ = True
__DEBUG__ = False
__TRACE__ = False
__KERNEL_DEBUG__ = False
-__PROFILE__ = True
+__PROFILE__ = False
__ENABLE_LONG_TESTS__ = False
diff --git a/hysop/backend/host/host_buffer.py b/hysop/backend/host/host_buffer.py
index c95a3b133..733c5f09b 100644
--- a/hysop/backend/host/host_buffer.py
+++ b/hysop/backend/host/host_buffer.py
@@ -8,7 +8,7 @@ class HostBuffer(np.ndarray, Buffer):
"""
Host buffer class.
"""
- __array_priority__ = +1.0
+ __array_priority__ = -1.0
def __new__(cls, size,
shape=None, dtype=np.uint8, order=None,
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index 5e8c81e90..95fff0992 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -17,6 +17,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
@debug
def __init__(self, **kwds):
super(PythonDirectionalAdvection, self).__init__(**kwds)
+ self.is_periodic = False
def get_work_properties(self):
requests = super(PythonDirectionalAdvection, self).get_work_properties()
@@ -55,11 +56,61 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
def setup(self, work):
super(PythonDirectionalAdvection, self).setup(work=work)
+ self.dpos = self.dposition.handle.view(type=npw.ndarray)
self.drtmp = work.get_buffer(self, 'rtmp', handle=True)
self.ditmp = work.get_buffer(self, 'itmp', handle=True)
self.dixtmp = work.get_buffer(self, 'ixtmp', handle=True)
if self.is_inplace:
self.dstmp = work.get_buffer(self, 'stmp', handle=True)
+
+ self._prepare_apply()
+
+ def _prepare_apply(self):
+ cr = self.compute_granularity
+ sinputs = self.dadvected_fields_in
+ soutputs = self.dadvected_fields_out
+
+ velo_mesh = self.dvelocity.mesh
+ velo_mesh_iterator = velo_mesh.build_compute_mesh_iterator(cr)
+ X0 = velo_mesh.local_compute_coords[-1]
+ dx = velo_mesh.space_step[-1]
+ inv_dx = (1.0/dx)
+ velo_compute_view = velo_mesh.local_compute_slices
+ self._velocity_mesh_attributes = (velo_mesh_iterator, dx, inv_dx, velo_compute_view, X0)
+
+ scalar_mesh = sinputs.values()[0].mesh
+ scalar_mesh_iterator = scalar_mesh.build_compute_mesh_iterator(cr)
+ sdx = scalar_mesh.space_step[-1]
+ inv_sdx = (1.0/sdx)
+ scalar_compute_view = scalar_mesh.local_compute_slices
+ self._scalar_mesh_attributes = (scalar_mesh_iterator, sdx, inv_sdx, scalar_compute_view)
+
+ remesh_ghosts = self.remesh_ghosts
+ outer_ghosts, inner_ghosts, proc_shape = {}, {}, {}
+ for sout in soutputs.values():
+ outer_ghosts[sout] = sout.get_outer_ghost_slices(remesh_ghosts)
+ inner_ghosts[sout] = sout.get_inner_ghost_slices(remesh_ghosts)
+ proc_shape[sout] = sout.mesh.proc_shape
+ self._outer_ghosts = outer_ghosts
+ self._inner_ghosts = inner_ghosts
+ self._proc_shape = proc_shape
+
+ in_compute_slices, out_compute_slices = {}, {}
+ in_ghosts, out_ghosts = {}, {}
+ in_shapes, out_shapes = {}, {}
+
+ for field in sinputs.keys():
+ Sin = sinputs[field]
+ Sout = soutputs[field]
+ in_compute_slices[Sin] = Sin.compute_slices
+ in_ghosts[Sin] = Sin.ghosts
+ in_shapes[Sin] = Sin.shape
+ out_compute_slices[Sout] = Sout.compute_slices
+ out_ghosts[Sout] = Sout.ghosts
+ out_shapes[Sout] = Sout.shape
+ self._inout_compute_slices = (in_compute_slices, out_compute_slices)
+ self._inout_ghosts = (in_ghosts, out_ghosts)
+ self._inout_shapes = (in_shapes, out_shapes)
@op_apply
def apply(self, simulation, **kwds):
@@ -111,22 +162,17 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
Vout[...] += Vin[Iy+(lidx,)]
def _compute_advection(self, dt):
- cr = self.compute_granularity
- P = self.dposition
- Vd = self.dvelocity[self.splitting_direction]
+ P = self.dpos
+ Vd = self.dvelocity(self.splitting_direction)
rtmp = self.drtmp
itmp = self.ditmp
- mesh = self.dvelocity.mesh
- X0 = mesh.local_compute_coords[-1]
- dx = mesh.space_step[-1]
- view = mesh.local_compute_slices
- dim = mesh.dim
-
- rk_scheme = self.time_integrator
+ cr = self.compute_granularity
+ is_periodic = self.is_periodic
+ rk_scheme = self.time_integrator
- is_periodic = False
- inv_dx = (1.0/dx)
+ (mesh_it,dx,inv_dx,view,X0) = self._velocity_mesh_attributes
+
_Vd = Vd[view]
Vd = Vd[view[:-1]+(slice(None),)]
@@ -140,7 +186,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
assert Vinf*dt*inv_dx <= cfl, msg
# fill position with first direction coordinates on the whole compute domain
- for idx,_,I,Ig in mesh.iter_compute_axes(cr):
+ for (idx,_,I,Ig) in mesh_it.iter_compute_mesh():
Pi = P[idx]
Vi = Vd[idx]
_Vi = _Vd[idx] # Vi without ghosts
@@ -202,31 +248,25 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
assert Pmax <= X0[-1] + Vmax*dt + 10*eps, msg
def _compute_remesh(self, dt):
- pos = self.dposition.handle
+ pos = self.dpos
sinputs = self.dadvected_fields_in
soutputs = self.dadvected_fields_out
is_inplace = self.is_inplace
- rtmp = self.drtmp
- itmp = self.ditmp
- stmp = self.dstmp
+
+ rtmp = self.drtmp
+ itmp = self.ditmp
+ stmp = self.dstmp
ixtmp = self.dixtmp
cr = self.compute_granularity
- mesh = sinputs.values()[0].mesh
- dim = mesh.dim
- dx = mesh.space_step[-1]
- view = mesh.local_compute_slices
+ (mesh_it, dx, inv_dx, compute_view) = self._scalar_mesh_attributes
scalar_advection_ghosts = self.scalar_advection_ghosts
remesh_ghosts = self.remesh_ghosts
remesh_kernel = self.remesh_kernel
P = 1 + remesh_kernel.n/2
- G = scalar_advection_ghosts + P - 1
- assert (G == remesh_ghosts)
-
- is_periodic = False
- inv_dx = (1.0/dx)
+ is_periodic = self.is_periodic
R0,R1 = rtmp[:2]
I0 = itmp[0]
@@ -236,14 +276,18 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
is_periodic = False
N = R0.shape[-1]
+
+ (in_compute_slices, out_compute_slices) = self._inout_compute_slices
+ (in_ghosts, out_ghosts) = self._inout_ghosts
+ (in_shapes, out_shapes) = self._inout_shapes
- for (idx,_,I,_) in mesh.iter_compute_axes(cr):
+ for (idx,_,I,_) in mesh_it.iter_compute_mesh():
Pi = pos[idx]
R0[...] = Pi
R0[...] *= inv_dx
I0[...] = npw.floor(R0).astype(I0.dtype)
R0[...] -= I0
- R0[...] *= -1 # we need -alpha
+ R0[...] *= -1 # we need -alpha for the left point
if __debug__:
Imin, Imax = I0.min(), I0.max()
@@ -260,22 +304,23 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
I0[...] += 1
R0[...] += 1
R1[...] = remesh_kernel.gamma(R0)
+
sid=0
for field in sinputs.keys():
Sin = sinputs[field]
Sout = soutputs[field]
- sin_view = tuple(idx[i] + Sin.ghosts[i] for i in xrange(cr))
- sin_view += Sin.compute_slices[cr:]
+ sin_view = tuple(idx[i] + in_ghosts[Sin][i] for i in xrange(cr))
+ sin_view += in_compute_slices[Sin][cr:]
- dG = Sout.ghosts[-1] - remesh_ghosts
- sout_view = tuple(idx[i] + Sout.ghosts[i] for i in xrange(cr))
- sout_view += Sout.compute_slices[cr:-1]
- sout_view += (slice(dG, Sout.shape[-1]-dG),)
+ dG = out_ghosts[Sout][-1] - remesh_ghosts
+ sout_view = tuple(idx[i] + out_ghosts[Sout][i] for i in xrange(cr))
+ sout_view += out_compute_slices[Sout][cr:-1]
+ sout_view += (slice(dG, out_shapes[Sout][-1]-dG),)
for k in xrange(field.nb_components):
- sin = sinputs[field].data[k].handle[sin_view]
- sout = soutputs[field].data[k].handle[sout_view]
+ sin = Sin(k)[sin_view]
+ sout = Sout(k)[sout_view]
Si = (S[sid] if is_inplace else sout)
if (q==-P+1):
Si[...] = 0.0
@@ -296,9 +341,9 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
# exchange and sum remeshed ghosts
for sout in soutputs.values():
buffers = sout.buffers
- outer_ghosts = sout.get_outer_ghost_slices(remesh_ghosts)
- inner_ghosts = sout.get_inner_ghost_slices(remesh_ghosts)
- nprocs = sout.mesh.proc_shape
+ outer_ghosts = self._outer_ghosts[sout]
+ inner_ghosts = self._inner_ghosts[sout]
+ nprocs = self._proc_shape[sout]
for buf in buffers:
lsrc, rsrc = outer_ghosts[-1]
@@ -312,7 +357,6 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
sout.exchange_ghosts()
-
@debug
def handle_method(self, method):
super(PythonDirectionalAdvection, self).handle_method(method)
diff --git a/hysop/backend/host/python/operator/transpose.py b/hysop/backend/host/python/operator/transpose.py
index 3c43b2a38..c5c53796b 100644
--- a/hysop/backend/host/python/operator/transpose.py
+++ b/hysop/backend/host/python/operator/transpose.py
@@ -25,13 +25,13 @@ class PythonTranspose(TransposeOperatorBase, HostOperator):
axes = self.axes
if self.is_inplace:
- din, dout, dtmp = self.din, self.dout, self.dtmp
+ din, dout, dtmp = self.din, self.dout, self.dtmp.handle
assert self.din.dfield is self.dout.dfield
for i in xrange(din.nb_components):
- dtmp[...] = np.transpose(din[i], axes=axes)
- dout[i][...] = dtmp
+ dtmp[...] = np.transpose(din(i), axes=axes)
+ dout(i)[...] = dtmp
else:
din, dout = self.din, self.dout
assert self.din.dfield is not self.dout.dfield
for i in xrange(din.nb_components):
- dout[i][...] = np.transpose(din[i], axes=axes)
+ dout(i)[...] = np.transpose(din(i), axes=axes)
diff --git a/hysop/core/graph/computational_graph.py b/hysop/core/graph/computational_graph.py
index ce06b335a..fd8999c31 100644
--- a/hysop/core/graph/computational_graph.py
+++ b/hysop/core/graph/computational_graph.py
@@ -48,6 +48,7 @@ class ComputationalGraph(ComputationalGraphNode):
self.nodes = []
self.graph = None
self.graph_built = False
+ self.graph_is_rendering = False
@graph_built
def _get_profiling_info(self):
@@ -339,6 +340,9 @@ class ComputationalGraph(ComputationalGraphNode):
win.connect("delete_event", Gtk.main_quit)
win.show_all()
Gtk.main()
+ self.graph_is_rendering = False
+
+ self.graph_is_rendering = True
from threading import Thread
display_thread = Thread(target=draw)
@@ -433,23 +437,25 @@ class ComputationalGraph(ComputationalGraphNode):
@debug
@ready
def apply(self, **kwds):
+ drawing = self.graph_is_rendering
reduced_graph = self.reduced_graph
operators = reduced_graph.vertex_properties['operators']
- active_ops = reduced_graph.vertex_properties['active_ops']
-
- old_color = None
- for vid in self.sorted_nodes:
- if old_color:
+ if drawing:
+ active_ops = reduced_graph.vertex_properties['active_ops']
+ old_color = None
+ for vid in self.sorted_nodes:
+ if old_color:
+ active_ops[vertex] = old_color
+ vertex = reduced_graph.vertex(vid)
+ old_color = active_ops[vertex]
+ active_ops[vertex] = 'red'
+ op = operators[vertex]
+ op.apply(**kwds)
active_ops[vertex] = old_color
- vertex = reduced_graph.vertex(vid)
- old_color = active_ops[vertex]
- active_ops[vertex] = 'red'
-
- op = operators[vertex]
- op_kwds = kwds.copy()
- op.apply(**op_kwds)
- active_ops[vertex] = old_color
+ else:
+ for op in self.nodes:
+ op.apply(**kwds)
@debug
@ready
diff --git a/hysop/core/graph/graph.py b/hysop/core/graph/graph.py
index 7f68ab59f..3c138f8d0 100644
--- a/hysop/core/graph/graph.py
+++ b/hysop/core/graph/graph.py
@@ -25,7 +25,7 @@ if __debug__:
def not_initialized(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _not_initialized(*args,**kargs):
return f(*args,**kargs)
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
@@ -33,12 +33,12 @@ if __debug__:
if self.initialized:
reason='this node has already been initialized.'
raise RuntimeError(msg.format(reason))
- return check
+ return _not_initialized
def initialized(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _initialized(*args,**kargs):
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
.format(self.__class__.__name__,f.__name__,self.name,'{}')
@@ -46,12 +46,12 @@ if __debug__:
reason='this node has not been initialized yet.'
raise RuntimeError(msg.format(reason))
return f(*args,**kargs)
- return check
+ return _initialized
def discretized(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _discretized(*args,**kargs):
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
.format(self.__class__.__name__,f.__name__,self.name,'{}')
@@ -59,12 +59,12 @@ if __debug__:
reason='this node has not been discretized yet.'
raise RuntimeError(msg.format(reason))
return f(*args,**kargs)
- return check
+ return _discretized
def ready(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _ready(*args,**kargs):
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
.format(self.__class__.__name__,f.__name__, self.name,'{}')
@@ -72,12 +72,12 @@ if __debug__:
reason='this node has not been set up.'
raise RuntimeError(msg.format(reason))
return f(*args,**kargs)
- return check
+ return _ready
def graph_built(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _graph_built(*args,**kargs):
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
.format(self.__class__.__name__,f.__name__,self.name,'{}')
@@ -85,12 +85,12 @@ if __debug__:
reason = 'the graph has not been built yet.'
raise RuntimeError(msg.format(reason))
return f(*args,**kargs)
- return check
+ return _graph_built
def generated(f):
assert callable(f)
@wraps(f)
- def check(*args,**kargs):
+ def _generated(*args,**kargs):
self = args[0]
msg = 'Cannot call {}.{}() on node \'{}\' because {}'\
.format(self.__class__.__name__,f.__name__,self.name,'{}')
@@ -98,7 +98,7 @@ if __debug__:
reason='this node has not been generated yet.'
raise RuntimeError(msg.format(reason))
return f(*args,**kargs)
- return check
+ return _generated
else: # not __debug__
# python optimized, no checks
diff --git a/hysop/core/memory/memory_request.py b/hysop/core/memory/memory_request.py
index 9112c12e9..3efc91348 100644
--- a/hysop/core/memory/memory_request.py
+++ b/hysop/core/memory/memory_request.py
@@ -327,7 +327,11 @@ class MultipleOperatorMemoryRequests(object):
raise ValueError(msg)
buffers = op_buffers[request_identifier]
if handle:
- buffers = tuple([b.handle for b in buffers])
+ from hysop.backend.host.host_buffer import HostBuffer
+ if isinstance(buffers[0].handle, HostBuffer):
+ buffers = tuple([b.handle.view(np.ndarray) for b in buffers])
+ else:
+ buffers = tuple([b.handle for b in buffers])
return buffers
def __str__(self):
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index e5296eb48..0bd1ac7d8 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -170,12 +170,15 @@ class CartesianDiscreteFieldView(DiscreteFieldView):
By default return all components.
"""
data = self.__get_data_views(*components)
- return tuple(d.handle for d in data)
+ if (self.backend.kind == Backend.HOST):
+ return tuple(d.handle.view(type=np.ndarray) for d in data)
+ else:
+ return tuple(d.handle for d in data)
def __getitem__(self, key):
"""
- Access to the content of the field.
- Returns ArrayBackend.Array component key of the field if key in an integer.
+ Access to the content of the field as with self.data.
+ Returns ArrayBackend.Array component(s) key of the field if key in an integer.
Else if key is a slice, returns a tuple of slice Arrays.
"""
if isinstance(key, slice):
@@ -183,10 +186,30 @@ class CartesianDiscreteFieldView(DiscreteFieldView):
return self._get_data(*range(*indices))
elif isinstance(key, (int,long)):
return self._get_data(key)[0]
+ elif isinstance(key, Ellipsis):
+ return self._get_data()
else:
msg='Unknown key type {} for CartesianDiscreteField.__getitem__().'
msg=msg.format(key.__class__)
raise TypeError(msg)
+
+ def __call__(self, key):
+ """
+ Access to the content of the field as with self.buffer.
+ Returns ArrayBackend.Array handle component(s) key of the field if key in an integer.
+ Else if key is a slice, returns a tuple of slice Arrays.
+ """
+ if isinstance(key, slice):
+ indices = key.indices(self.nb_components)
+ return self._get_buffers(*range(*indices))
+ elif isinstance(key, (int,long)):
+ return self._get_buffers(key)[0]
+ elif isinstance(key, Ellipsis):
+ return self._get_buffers()
+ else:
+ msg='Unknown key type {} for CartesianDiscreteField.__call__().'
+ msg=msg.format(key.__class__)
+ raise TypeError(msg)
def _get_axes(self):
"""Return the permutation scheme in numpy notations."""
diff --git a/hysop/mesh/cartesian_mesh.py b/hysop/mesh/cartesian_mesh.py
index dbda4d668..dd78b0ea1 100644
--- a/hysop/mesh/cartesian_mesh.py
+++ b/hysop/mesh/cartesian_mesh.py
@@ -352,20 +352,41 @@ class CartesianMeshView(MeshView):
for idx in np.ndindex(*iter_shape):
yield (idx, I)
- def iter_compute_axes(self, axes=0):
+
+ def iter_compute_mesh(self, axes=0):
"""
Return an iterator over a compute mesh array that returns indices.
Iterates 'axes' axes at a time (0=whole array view,
1=one axe at a time, max value is self.dim).
"""
+ return self.build_compute_mesh_iterator(axes=axes).iter_compute_mesh()
+
+ def build_compute_mesh_iterator(self, axes):
+ """
+ Return an object that can iterate over the local
+ computational mesh multiple times.
+ Usefull for optimizing operator apply methods.
+ See self.iter_compute_mesh().
+ """
assert 0 <= axes <= self.dim, 'bad value for axes'
- ghosts = self.ghosts
- iter_shape = tuple(self.compute_resolution[:axes])
- gI = np.ix_(*self.local_compute_indices[axes:])
- for idx in np.ndindex(*iter_shape):
- gidx = tuple(idx[i]+ghosts[i] for i in xrange(axes))
- I = tuple(gI[i]-ghosts[j] for i,j in enumerate(xrange(axes, self.dim)))
- yield (idx, gidx, I, gI)
+
+ class __CartesianMeshComputeAxeIterator(object):
+ def __init__(self, dim, ghosts, compute_resolution, local_compute_indices):
+ iter_shape = tuple(compute_resolution[:axes])
+ gI = np.ix_(*local_compute_indices[axes:])
+ self._data = (dim, ghosts, iter_shape, gI)
+
+ def iter_compute_mesh(self):
+ (dim, ghosts, iter_shape, gI) = self._data
+ for idx in np.ndindex(*iter_shape):
+ gidx = tuple(idx[i]+ghosts[i] for i in xrange(axes))
+ I = tuple(gI[i]-ghosts[j] for i,j in enumerate(xrange(axes, dim)))
+ yield (idx, gidx, I, gI)
+
+ return __CartesianMeshComputeAxeIterator(dim=self.dim,
+ ghosts=self.ghosts,
+ compute_resolution=self.compute_resolution,
+ local_compute_indices=self.local_compute_indices)
def is_inside_local_domain(self, point):
"""
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 17b7c2bf2..c2e436c06 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -100,32 +100,39 @@ class CflTimestepCriteria(TimestepCriteria):
super(CflTimestepCriteria,self).__init__(name=name, variables={velocity:None} )
self.cfl = cfl
self.velocity = velocity
+ self._compute_criteria_impl = None
def discretize(self):
super(CflTimestepCriteria, self).discretize()
- self.dvelocity = self.input_discrete_fields[self.velocity]
-
- def compute_criteria(self, **kwds):
+ dv = self.input_discrete_fields[self.velocity]
+ self.dvelocity = dv
+
from hysop.fields.cartesian_discrete_field import CartesianDiscreteFieldView
-
- dv = self.dvelocity
- cfl = self.cfl
- dt = +npw.inf
-
if isinstance(dv, CartesianDiscreteFieldView):
- dx = dv.dfield.mesh.space_step
- 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))
- if Vi_inf == 0:
- dti = +npw.inf
- else:
- dti = (dx[i]*cfl) / Vi_inf
- dt = min(dt, dti)
+ self._compute_criteria_impl = self._compute_criteria_cartesian
+ self._dx = dv.dfield.mesh.space_step
else:
msg='CFL criteria has not been implemented for {} discrete fields.'
msg=msg.format(type(self.dvelocity))
raise RuntimeError(msg)
+
+ def compute_criteria(self, **kwds):
+ return self._compute_criteria_impl(**kwds)
+
+ def _compute_criteria_cartesian(self, **kwds):
+ dv = self.dvelocity
+ cfl = self.cfl
+ dt = +npw.inf
+ dx = self._dx
+
+ 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))
+ if Vi_inf == 0:
+ dti = +npw.inf
+ else:
+ dti = (dx[i]*cfl) / Vi_inf
+ dt = min(dt, dti)
assert dt>0.0
return dt
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index f40541064..eb3f41159 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -15,6 +15,7 @@ from hysop.tools.decorators import debug
from hysop.tools.types import check_instance
from hysop.tools.numpywrappers import npw
from hysop.tools.io_utils import IO, IOParams, XMF
+from hysop.core.graph.graph import op_apply
from hysop.core.graph.computational_graph import ComputationalGraphOperator
from hysop.fields.continuous_field import Field
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
@@ -127,6 +128,7 @@ class HDF_IO(ComputationalGraphOperator):
refmesh = self.subset.mesh[self.topology]
else:
refmesh = self.topology.mesh
+ self.refmesh = refmesh
# Global resolution for hdf5 output (warning : this must
# be the whole domain resolution, not the subset resolution)
@@ -211,10 +213,8 @@ class HDF_Writer(HDF_IO):
if xmfalways:
self.step = self._step_HDF5_XMF
- self.finalize = lambda: 1
else:
self.step = self._step_HDF5
- self.finalize = self.createXMFFile
self._xdmf_data_files = []
# filename = prefix_N, N = counter value
self._get_filename = self._input_fname
@@ -223,6 +223,17 @@ class HDF_Writer(HDF_IO):
# time where xdmf has been written, to raise an exception
# if that happens.
self._last_written_time = None
+ self._xmf_file = None
+
+ def setup(self, **kwds):
+ super(HDF_Writer,self).setup(**kwds)
+ self._setup_grid_template()
+ self.createXMFFile()
+
+ def finalize(self):
+ self.updateXMFFile()
+ if self._xmf_file:
+ self._xmf_file.close()
def _input_fname(self, i):
"""Set output file name for current iteration"""
@@ -230,7 +241,7 @@ class HDF_Writer(HDF_IO):
assert i >= 0, msg
return self.io_params.filename + "_{0:05d}".format(i) + '.h5'
- @debug
+ @op_apply
def apply(self, simulation=None, **kwds):
if simulation is None:
raise ValueError("Missing simulation value for monitoring.")
@@ -241,11 +252,41 @@ class HDF_Writer(HDF_IO):
if self.kill:
sys.exit(1)
- def createXMFFile(self):
- """Create and fill the xdmf file
- """
+ def _setup_grid_template(self):
+ topo = self.topology
+ dim = topo.domain.dim
+ dx = list(topo.mesh.space_step)
+ mesh = self.refmesh
+ subset = self.subset
+ if subset is not None:
+ res = list(mesh[topo].grid_resolution)
+ orig = list(subset.real_orig[topo])
+ else:
+ res = list(mesh.grid_resolution)
+ orig = list(topo.domain.origin)
+ resolution = [1,]*3
+ origin = [0.0,]*3
+ step = [0.0,]*3
+
+ idim = 3-dim
+ resolution[idim:] = res
+ origin[idim:] = orig
+ step[idim:] = dx
+
+ write_resolution = tuple(resolution)
+ write_origin = tuple(origin)
+ write_step = tuple(step)
+
+ ds_names = self.dataset.keys()
+ grid_attributes = XMF.prepare_grid_attributes(
+ ds_names,
+ resolution, origin, step)
+ self.grid_attributes_template = grid_attributes
+
- if self.mpi_params.rank == self.io_params.io_leader:
+ def createXMFFile(self):
+ """Create and fill the header of the xdmf file."""
+ if (self.mpi_params.rank == self.io_params.io_leader):
f = open(self.io_params.filename + '.xmf', 'w')
f.write("<?xml version=\"1.0\" ?>\n")
f.write("<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\">\n")
@@ -253,15 +294,29 @@ class HDF_Writer(HDF_IO):
f.write(" <Domain>\n")
f.write(" <Grid Name=\"CellTime\" GridType=\"Collection\" ")
f.write("CollectionType=\"Temporal\">\n")
- ds_names = self.dataset.keys()
- for i, t in self._xdmf_data_files:
- f.write(XMF.write_grid_attributes(
- self.topology, ds_names, i, t, self._get_filename(i),
- self.subset))
- f.write(" </Grid>\n")
- f.write(" </Domain>\n")
- f.write("</Xdmf>\n")
- f.close()
+ self._last_xmf_pos = f.tell()
+ f.flush()
+ self._xmf_file = f
+
+ def updateXMFFile(self):
+ """Update xdmf file."""
+ if (self.mpi_params.rank == self.io_params.io_leader):
+ f = self._xmf_file
+ lastp = self._last_xmf_pos
+ assert (f is not None)
+ assert (lastp is not None)
+ for (i, t) in self._xdmf_data_files:
+ filename = self._get_filename(i).split('/')[-1]
+ grid_attrs = self.grid_attributes_template.format(
+ niteration=i, time=t, filename=filename)
+ f.seek(lastp)
+ f.write(grid_attrs)
+ self._last_xmf_pos = f.tell()
+ f.write(" </Grid>\n")
+ f.write(" </Domain>\n")
+ f.write("</Xdmf>\n")
+ f.flush()
+ self._xdmf_data_files = []
def _step_HDF5(self, simu):
"""Write an h5 file with data on each mpi process.
@@ -299,8 +354,7 @@ class HDF_Writer(HDF_IO):
def _step_HDF5_XMF(self, simu):
self._step_HDF5(simu)
- self.createXMFFile()
-
+ self.updateXMFFile()
class HDF_Reader(HDF_IO):
"""
@@ -335,8 +389,8 @@ class HDF_Reader(HDF_IO):
self.io_params.filename + "_{0:05d}".format(i) + '.h5'
else:
self._get_filename = lambda i=None: self.io_params.filename
-
- @debug
+
+ @op_apply
def apply(self, simulation=None, **kwds):
# Read HDF file
self.open_hdf(count=self.restart, mode='r')
diff --git a/hysop/parameters/numpy_parameter.py b/hysop/parameters/numpy_parameter.py
index e1dc61c1a..b0d73aaf0 100644
--- a/hysop/parameters/numpy_parameter.py
+++ b/hysop/parameters/numpy_parameter.py
@@ -153,11 +153,11 @@ class NumpyParameter(Parameter):
assert (value is not None) and isinstance(value, np.ndarray)
if (value.shape != self.shape):
msg='Parameter shape mismatch, expected {} but got {}.'
- msg=msg.format(self.shape, other.shape)
+ msg=msg.format(self.shape, value.shape)
raise ValueError(msg)
if (value.dtype != self.dtype):
msg='Parameter dtype mismatch, expected {} but got {}.'
- msg=msg.format(self.dtype, other.dtype)
+ msg=msg.format(self.dtype, value.dtype)
raise ValueError(msg)
if (self.min_value is not None) or \
(self.max_value is not None) or \
diff --git a/hysop/parameters/scalar_parameter.py b/hysop/parameters/scalar_parameter.py
index a82f95ac7..bb3c3146f 100644
--- a/hysop/parameters/scalar_parameter.py
+++ b/hysop/parameters/scalar_parameter.py
@@ -26,7 +26,8 @@ class ScalarParameter(NumpyParameter):
"""
if not isinstance(value, np.ndarray):
assert np.isscalar(value), 'value is neither a np.ndarray, nor a scalar.'
- value = np.asarray([value], dtype=self.dtype).reshape((1,))
+ value = np.asarray([value], dtype=self.dtype)
+ value = value.reshape((1,))
super(ScalarParameter,self)._set_value_impl(value)
def __call__(self):
diff --git a/hysop/problem.py b/hysop/problem.py
index 2daf17286..43a8dcad0 100644
--- a/hysop/problem.py
+++ b/hysop/problem.py
@@ -26,7 +26,7 @@ class Problem(ComputationalGraph):
self.setup(work)
@debug
- def solve(self, simu, max_iter=None, report_freq=10, **kargs):
+ def solve(self, simu, max_iter=None, report_freq=20, **kargs):
vprint('\nSolving problem...')
simu.initialize()
while not simu.is_over:
diff --git a/hysop/tools/decorators.py b/hysop/tools/decorators.py
index 255a621c1..1ae124dcf 100644
--- a/hysop/tools/decorators.py
+++ b/hysop/tools/decorators.py
@@ -124,13 +124,13 @@ def profile(f):
if __PROFILE__:
from hysop.core.mpi import Wtime
@wraps(f)
- def decorator(*args, **kwargs):
+ def _profile(*args, **kwargs):
"""args[0] contains the object"""
t0 = Wtime()
res = f(*args, **kwargs)
args[0]._profiler[f.func_name] += Wtime() - t0
return res
- return decorator
+ return _profile
else:
return f
diff --git a/hysop/tools/io_utils.py b/hysop/tools/io_utils.py
index dcdc35af4..3211a174b 100755
--- a/hysop/tools/io_utils.py
+++ b/hysop/tools/io_utils.py
@@ -419,23 +419,18 @@ class XMF(object):
return buff.replace(']', '').replace('(', '').replace(')', '')
@staticmethod
- def write_grid_attributes(topo, dataset_names, ite,
- time, filename, subset=None):
+ def prepare_grid_attributes(dataset_names,
+ resolution, origin, step):
"""
- Write XDMF header into a file
+ Prepare XDMF header as a string.
Parameters
-----------
- topo : :class:`hysop.topology.topology.CartesianTopology`
- used as reference to define local and global meshes in xdmf file.
dataset_names : list
all datasets names
- ite : int
- iteration number
- time : double
- current time
- filename : string
- name of the hdf file which contains datas for the current process.
+ resolution: 3d tuple
+ origin: 3d tuple
+ step: 3d tuple
subset : :class:`hysop.domain.subsets.Subset`, optional
to define a grid only on this subset.
If None, grid on the whole domain (from topo)
@@ -443,38 +438,22 @@ class XMF(object):
Returns:
--------
string
- the xml-like header.
+ the xml-like header formattable with the following keywords:
+ niteration : iteration number
+ time: time in seconds
+ filename: target file name
"""
# The header (xml-like), saved in a string.
# always use a 3D mesh because paraview does not like 2D meshes (uses axe (Y,Z) instead of (X,Y)).
xml_grid = ""
- dim = topo.domain.dim
topo_type = "3DCORECTMesh"
geo_type = "ORIGIN_DXDYDZ"
- xml_grid += " <Grid Name=\"Iteration {0:03d}\"".format(ite)
+ xml_grid += " <Grid Name=\"Iteration {}\"".format('{niteration:03d}')
xml_grid += " GridType=\"Uniform\">\n"
- xml_grid += " <Time Value=\"{0}\" />\n".format(time)
+ xml_grid += " <Time Value=\"{}\" />\n".format('{time}')
xml_grid += " <Topology TopologyType=\"" + str(topo_type) + "\""
xml_grid += " NumberOfElements=\""
-
- # Check substet to find the required grid resolution
- dx = list(topo.mesh.space_step)
- if subset is not None:
- res = list(subset.mesh[topo].grid_resolution)
- orig = list(subset.real_orig[topo])
- else:
- res = list(topo.mesh.grid_resolution)
- orig = list(topo.domain.origin)
- resolution = [1,]*3
- origin = [0.0,]*3
- step = [0.0,]*3
-
- idim = 3-dim
- resolution[idim:] = res
- origin[idim:] = orig
- step[idim:] = dx
-
xml_grid += XMF._list_format(resolution) + " \"/>\n"
xml_grid += " <Geometry GeometryType=\"" + geo_type + "\">\n"
xml_grid += " <DataItem Dimensions=\"" + str(3) + " \""
@@ -495,7 +474,7 @@ class XMF(object):
xml_grid += XMF._list_format(resolution) + " \""
xml_grid += " NumberType=\"Float\" Precision=\"8\" Format=\"HDF\""
xml_grid += " Compression=\"Raw\">\n" #
- xml_grid += " " + filename.split('/')[-1]
+ xml_grid += " {filename}"
xml_grid += ":/" + name
xml_grid += "\n </DataItem>\n"
xml_grid += " </Attribute>\n"
--
GitLab