From cfbdf41da8106d13c299d71a5e306599a9c2be37 Mon Sep 17 00:00:00 2001
From: Keck Jean-Baptiste <jean-baptiste.keck@imag.fr>
Date: Fri, 22 Sep 2017 23:09:50 +0200
Subject: [PATCH] fixed some bugs, added dfield.clone() method, directional
 advection testing

---
 hysop/__init__.py                             |   2 +-
 .../operator/directional/advection_dir.py     |   6 +-
 hysop/core/graph/graph_builder.py             |   6 +-
 hysop/core/memory/memory_request.py           |   3 +-
 hysop/deps.py                                 |   2 +-
 hysop/fields/cartesian_discrete_field.py      |  69 +++++++--
 hysop/fields/discrete_field.py                |  30 +++-
 .../tests/test_directional_advection.py       | 140 +++++++++++-------
 hysop/parameters/numpy_parameter.py           |  20 +--
 hysop/parameters/parameter.py                 |  43 +++---
 hysop/simulation.py                           |  12 +-
 hysop/tools/transposition_states.py           |   5 +-
 hysop/topology/cartesian_topology.py          |   2 +
 13 files changed, 221 insertions(+), 119 deletions(-)

diff --git a/hysop/__init__.py b/hysop/__init__.py
index 3afa33583..ed9bf2cc3 100644
--- a/hysop/__init__.py
+++ b/hysop/__init__.py
@@ -16,7 +16,7 @@ __FFTW_ENABLED__   = "ON" is "ON"
 __SCALES_ENABLED__ = "ON" is "ON"
 __OPTIMIZE__       = not __debug__
 
-__VERBOSE__        = False
+__VERBOSE__        = True
 __DEBUG__          = False
 __TRACE__          = False
 __KERNEL_DEBUG__   = False
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index c4db1f92c..ebe356aa6 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -47,8 +47,8 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
         if self.is_inplace:
             rg = self.remesh_ghosts
             sshape = shape[:-1] + (shape[-1]+2*rg,)
-            request = MemoryRequest.empty_like(a=V[-1], shape=sshape, 
-                    nb_components=nscalars)
+            request = MemoryRequest.empty_like(a=V, shape=sshape, 
+                                                nb_components=nscalars)
             requests.push_mem_request('stmp', request)
         self.nscalars = nscalars
 
@@ -299,7 +299,6 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
         last_scalar_axe_mesh_indices = self._last_scalar_axe_mesh_indices
         
         for (idx,_,I,_) in mesh_it.iter_compute_mesh():
-            print 'LOLOLOLO', idx
             Pi = pos[idx]
             R0[...]  = Pi
             R0[...] *= inv_dx
@@ -362,6 +361,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
                         if (q==P) and is_inplace:
                             sout[...] = Si
                         sid+=1
+
         # exchange and sum remeshed ghosts
         ghost_exchangers = self._out_ghost_exchangers
         for sout in dsoutputs.values():
diff --git a/hysop/core/graph/graph_builder.py b/hysop/core/graph/graph_builder.py
index dbe41748f..484e24615 100644
--- a/hysop/core/graph/graph_builder.py
+++ b/hysop/core/graph/graph_builder.py
@@ -933,11 +933,12 @@ class GraphBuilder(object):
                         if target_dfield_requirements:
                             tstates = target_dfield_requirements.transposition_states
                             default_tstate = TranspositionState[dim].default()
+                            print 'before:', istate.transposition_state
                             if (tstates is None) or (default_tstate in tstates):
                                 istate.transposition_state = default_tstate
                             else:
-                                print tstates
                                 istate.transposition_state = tstates[0]
+                            print 'after:', istate.transposition_state
 
                             basis = target_dfield_requirements.basis
                             default_basis = (Basis.default,)*dim
@@ -946,6 +947,9 @@ class GraphBuilder(object):
                             else:
                                 istate.basis = basis[0]
                             dprint('       >Initial state set to {}'.format(istate))
+                        print 'TARGET TSTATES: ', tstates 
+                        print 'DEFAULT TSTATE: ', default_tstate
+                        print 'ISTATE:\n', istate 
                     else:
                         istate = dtopology_states.setdefault(target_topo, 
                                                 CartesianTopologyState(dim))
diff --git a/hysop/core/memory/memory_request.py b/hysop/core/memory/memory_request.py
index 2facdf55d..d7f0c202e 100644
--- a/hysop/core/memory/memory_request.py
+++ b/hysop/core/memory/memory_request.py
@@ -345,7 +345,8 @@ class MultipleOperatorMemoryRequests(object):
                 precision = ''
             s+= ' Backend {}{}:'.format(kind, precision)
             total=0
-            for op,op_requests in backend_requests.iteritems():
+            for op in sorted(backend_requests.keys(), key=lambda op: op.name):
+                op_requests = backend_requests[op]
                 s+= '\n  Operator {} ({})'.format(op.name, op.__class__.__name__)
                 local_total=0
                 for req in op_requests:
diff --git a/hysop/deps.py b/hysop/deps.py
index bc5468d2f..16e26f6da 100644
--- a/hysop/deps.py
+++ b/hysop/deps.py
@@ -22,7 +22,7 @@ except ImportError as e:
 
 import sys, os, subprocess, platform, warnings, traceback
 import resource, psutil, tempfile, cpuinfo
-import inspect, functools, operator
+import inspect, functools, operator, random
 import hashlib, gzip, copy, types, string
 import math, re, contextlib
 import six, itertools
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index add6b4226..96cd75b18 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -9,8 +9,8 @@ Documentation and examples can be found in :ref:`fields`.
 from hysop import vprint, dprint
 from hysop.deps import np
 from hysop.constants import Backend
-from hysop.tools.decorators import debug
-from hysop.tools.types import check_instance, to_tuple
+from hysop.tools.decorators import debug, static_vars
+from hysop.tools.types import check_instance, to_tuple, first_not_None
 from hysop.tools.numpywrappers import npw
 from hysop.tools.misc import prod
 from hysop.fields.discrete_field import DiscreteField, DiscreteFieldView
@@ -286,20 +286,36 @@ CartesianDiscreteFieldView (id={}, tag={})
         self.shape, self.dtype, self.initial_values,
         self.topology.tag, self.topology_state.short_description())
         return s
-
-    def fill(self, value, **kwds):
-        """Set all components to specified value."""
-        for d in xrange(self.nb_components):
-            self.data[d].fill(value=value, **kwds)
     
-    def randomize(self, **kwds):
-        """Initialize a the with random values."""
-        for d in xrange(self.nb_components):
-            self.backend.rand(out=self.data[d], **kwds)
     
+    def clone(self, name=None, tstate=None):
+        """
+        Create a new temporary DiscreteField and allocate it 
+        like the current object, possibly on a different backend. 
+
+        This should only be used for debugging and testing purpose.
+        The generated discrete field is not registered to the continuous
+        field.
+        """
+        
+        default_name='{}_clone_{}'.format(self.name, self._dfield._clone_id)
+        self._dfield._clone_id += 1
+
+        tstate = first_not_None(tstate, self._topology_state)
+        name = first_not_None(name, default_name)
+        
+        dfield = CartesianDiscreteField(name=name, 
+                field=self._dfield._field, 
+                topology=self._dfield._topology, 
+                init_topology_state=tstate,
+                is_cloned=True)
+
+        dfield.copy(self)
+        return dfield
 
     def copy(self, field_in, **kwds):
-        """Initialize a field on topo with values from another field.
+        """
+        Initialize a field on topo with values from another field.
 
         Parameters
         -----------
@@ -308,6 +324,16 @@ CartesianDiscreteFieldView (id={}, tag={})
         """
         for d in xrange(self.nb_components):
             self.backend.memcpy(dst=self.data[d], src=field_in[d], **kwds)
+    
+    def fill(self, value, **kwds):
+        """Set all components to specified value."""
+        for d in xrange(self.nb_components):
+            self.data[d].fill(value=value, **kwds)
+    
+    def randomize(self, **kwds):
+        """Initialize a the with random values."""
+        for d in xrange(self.nb_components):
+            self.backend.rand(out=self.data[d], **kwds)
    
     def initialize(self, formula, vectorize=False, **kwds):
         """
@@ -366,7 +392,7 @@ CartesianDiscreteFieldView (id={}, tag={})
             vprint(msg)
             self.data = data
 
-    def norm(self, p):
+    def norm(self, p=2, data=None):
         """
         Compute the euclidian norm of the discrete field
         p-norm = (sum(data[d]**p)**1/p for d = 1..dim
@@ -375,13 +401,26 @@ CartesianDiscreteFieldView (id={}, tag={})
         result  = npw.zeros((self.nb_components))
         gresult = npw.zeros((self.nb_components))
         slices = self.compute_slices
-        data = self.data
+        data = first_not_None(data, self.data)
         for d in range(self.nb_components):
             tmp = abs(data[d][slices])**p
-            result[d] = self.array_backend.sum(tmp).get()
+            result[d] = tmp.sum().get()
         self.topology.cart_comm.Allreduce(result, gresult)
         return gresult ** (1.0/p)
 
+    def distance(self, other, p=2):
+        """
+        Compute the distance between two discrete fields
+        p-distance = (sum(|data[d]|**p)**1/p for d = 1..dim
+          where data[d] = |self.data[d] - other.data[d]|
+        sum on all grid points excluding ghosts.
+        """
+        check_instance(other, CartesianDiscreteField)
+        assert self.nb_components == other.nb_components
+        data = tuple(self.data[i] - other.data[i] for i in xrange(self.nb_components))
+        return self.norm(p=p, data=data)
+
+
     def print_with_ghosts(self, component=0, compute=None, 
             inner_ghosts=None, outer_ghosts='X', **print_opts):
         """
diff --git a/hysop/fields/discrete_field.py b/hysop/fields/discrete_field.py
index 286bc3e04..725c483cf 100644
--- a/hysop/fields/discrete_field.py
+++ b/hysop/fields/discrete_field.py
@@ -114,7 +114,19 @@ class DiscreteFieldView(TaggedObjectView, VariableTag):
         pass
     @abstractmethod
     def copy(self, from_dfield, **kwds):
-        """Copy to this discrete field from another one."""
+        """Fill this discrete field with values from another one."""
+        pass
+
+    @abstractmethod
+    def clone(self, tstate=None):
+        """
+        Create a new temporary DiscreteField and allocate it 
+        like the current object, with possibly a different topology state.
+
+        This should only be used for debugging and testing purpose.
+        The generated discrete field is not registered to the continuous
+        field.
+        """
         pass
 
     @abstractmethod
@@ -188,7 +200,7 @@ class DiscreteField(TaggedObject):
     __metaclass__ = ABCMeta
     
     @debug
-    def __init__(self, field, topology, name=None, **kwds):
+    def __init__(self, field, topology, name=None, is_cloned=False, **kwds):
         """
         Creates a discrete field for a given continuous field and topology.
 
@@ -213,16 +225,18 @@ class DiscreteField(TaggedObject):
             msg=msg.format(field.domain.tag, topology.domain.tag)
             raise ValueError(msg)
         
-        if topology in field.discrete_fields:
-            msg='Field {} was already discretized on topology {}.'.format(
-                    field.name, topology.tag)
-            raise RuntimeError(msg)
-
-        field.discrete_fields[topology] = self
+        if not is_cloned:
+            if topology in field.discrete_fields:
+                msg='Field {} was already discretized on topology {}.'.format(
+                        field.name, topology.tag)
+                raise RuntimeError(msg)
+        
+            field.discrete_fields[topology] = self
 
         self._name  = name or '{}{}'.format(field.name, topology.id)
         self._topology = topology
         self._field = field
+        self._clone_id = 0
     
     @abstractmethod
     def view(self, topology_state):
diff --git a/hysop/operator/tests/test_directional_advection.py b/hysop/operator/tests/test_directional_advection.py
index 12dbe13e4..bc4487f8e 100644
--- a/hysop/operator/tests/test_directional_advection.py
+++ b/hysop/operator/tests/test_directional_advection.py
@@ -1,6 +1,5 @@
 
-import random
-from hysop.deps import np, it
+from hysop.deps import np, it, sys, random
 from hysop.testsenv import __ENABLE_LONG_TESTS__, __HAS_OPENCL_BACKEND__
 from hysop.testsenv import opencl_failed, iter_clenv
 from hysop.tools.contexts import printoptions
@@ -10,6 +9,7 @@ from hysop.tools.io_utils import IO
 from hysop.operator.directional.advection_dir import DirectionalAdvection
 
 from hysop import Field, Box, Simulation
+from hysop.methods import Remesh, TimeIntegrator
 from hysop.constants import Implementation
 from hysop.numerics.splitting.strang import StrangSplitting, StrangOrder
 from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK3, RK4
@@ -28,7 +28,7 @@ class TestDirectionalAdvectionOperator(object):
             cls.size_max = 8
         else:
             cls.size_min = 8
-            cls.size_max = 32
+            cls.size_max = 16
         
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
@@ -42,23 +42,26 @@ class TestDirectionalAdvectionOperator(object):
             size_min=None, size_max=None):
         enable_extra_tests = self.enable_extra_tests
         assert dim > 0
-
+        
+        # periodic boundaries removes one computational point
+        # so we add one here.
         size_min = first_not_None(size_min, self.size_min) + 1
         size_max = first_not_None(size_max, self.size_max) + 1
         
         shape = tuple(np.random.randint(low=size_min, high=size_max, size=dim).tolist())
-        nb_components = 2
+        nb_scalar_components = 2
         
         domain = Box(length=(2*np.pi,)*dim)
-        for dtype in (np.float32, np.float64):
-            for velocity_cfl in (0.55, 1.75):
-                Vin  = Field(domain=domain, name='Vin', dtype=dtype,  
-                        nb_components=nb_components, register_object=False)
-                Sin  = Field(domain=domain, name='Sin', dtype=dtype,  
-                        nb_components=nb_components, register_object=False)
-                Sout = Field(domain=domain, name='Sout', dtype=dtype, 
-                        nb_components=nb_components, register_object=False)
+        for dtype in (np.float32, np.float64, np.float):
+            
+            Vin  = Field(domain=domain, name='Vin', dtype=dtype,  
+                    nb_components=dim, register_object=False)
+            Sin  = Field(domain=domain, name='Sin', dtype=dtype,  
+                    nb_components=nb_scalar_components, register_object=False)
+            Sout = Field(domain=domain, name='Sout', dtype=dtype, 
+                    nb_components=nb_scalar_components, register_object=False)
 
+            for velocity_cfl in (0.55, 1.75):
                 self._test_one(shape=shape, dim=dim, dtype=dtype,
                         is_inplace=is_inplace, domain=domain, 
                         Vin=Vin, Sin=Sin, Sout=Sout, velocity_cfl=velocity_cfl)
@@ -73,11 +76,12 @@ class TestDirectionalAdvectionOperator(object):
                 d[...] = 0.0
 
     @classmethod
-    def __scalar_init(cls, data, coords):
+    def __scalar_init(cls, data, coords, offsets=None):
+        offsets = first_not_None(offsets, (0.0,)*len(coords))
         for i,d in enumerate(data):
             d[...] = 1.0/(i+1)
-            for coord in coords:
-                d[...] *= np.cos(coord)
+            for coord, offset in zip(coords, offsets):
+                d[...] *= np.cos(coord+offset)
 
 
     def _test_one(self, shape, dim,
@@ -96,14 +100,15 @@ class TestDirectionalAdvectionOperator(object):
             variables = { vin: shape, sin: shape, sout: shape }
 
         implementations = DirectionalAdvection.implementations()
+        method = {TimeIntegrator: Euler, Remesh: Remesh.L4_4}
        
         def iter_impl(impl):
             base_kwds = dict(velocity=vin, advected_fields=sin, advected_fields_out=sout, 
                     velocity_cfl=velocity_cfl, variables=variables, implementation=impl, 
-                    name='test_advec_{}'.format(str(impl).lower()))
+                    method=method, name='test_advec_{}'.format(str(impl).lower()))
             
             if impl is Implementation.PYTHON:
-                yield 'default implementation: ', DirectionalAdvection(**base_kwds)
+                yield 'default', DirectionalAdvection(**base_kwds)
             elif impl is Implementation.OPENCL_CODEGEN:
                 # for cl_env in iter_clenv():
                     # msg='platform {}, device {}: '.format(cl_env.platform.name, 
@@ -115,39 +120,60 @@ class TestDirectionalAdvectionOperator(object):
                 raise NotImplementedError(msg)
         
         # Compare to other implementations
-        advec_axes = set(tuple((x,) for x in xrange(dim)) + (tuple(xrange(dim)),))
+        advec_axes = ((),)
+        advec_axes += tuple((x,) for x in xrange(dim))
+        if dim>1:
+            advec_axes += (tuple(xrange(dim)),)
         for impl in implementations:
             print ' >Implementation {}:'.format(impl)
             for sop,op in iter_impl(impl):
-                print '  *{}:'.format(sop)
+                print '  *{}: '.format(sop)
                 for axes in advec_axes:
-                    print '.'
-                    split = StrangSplitting(splitting_dim=dim, 
-                                            order=StrangOrder.STRANG_SECOND_ORDER)
-                    split.push_operators(op)
-                    split.build()
-
-                    dvin  = split.input_discrete_fields[vin]
-                    dsin  = split.input_discrete_fields[sin]
-                    dsout = split.output_discrete_fields[sout]
-                    
-                    dvin.initialize(self.__velocity_init, axes=axes)
-                    dsin.initialize(self.__scalar_init)
-
-                    dt = (0.95 * velocity_cfl) / (max(shape[i] for i in axes)-1)
-                    napplies = 5
-                    simu = Simulation(start=0.0, end=5*dt, dt0=dt, ensure_unique=False)
-                    simu.initialize()
-                    while not simu.is_over:
-                        split.apply(simulation=simu)
-                        simu.advance()
-                    
-                    out = tuple( data.get().handle for data in dsout.data )
-                    # dref  = CartesianDiscreteField.dfield_like(dsout, backend=Backend)
-                    # dref.initialize(__scalar_init, offset=dt*)
-
-                    simu.finalize()
-                    split.finalize()
+                    Vi = np.asarray([+1.0 if (i in axes) else +0.0 for i in xrange(dim)], dtype=dtype)
+                    #print 'Vi = {}'.format(Vi)
+                    #sys.stdout.write('.')
+                    #sys.stdout.flush()
+                    try:
+                        split = StrangSplitting(splitting_dim=dim, 
+                                                order=StrangOrder.STRANG_SECOND_ORDER)
+                        split.push_operators(op)
+                        split.build()
+                        split.display()
+                        sys.exit(1)
+
+                        dvin  = split.input_discrete_fields[vin]
+                        dsin  = split.input_discrete_fields[sin]
+                        dsout = split.output_discrete_fields[sout]
+                        dsref = dsout.clone()
+                        view  = dsout.compute_slices
+                        
+                        dvin.initialize(self.__velocity_init, axes=axes)
+                        dsin.initialize(self.__scalar_init)
+                        
+                        # Use optimal timestep, ||Vi||_inf is 1 on a per-axis basis
+                        dt = (0.95 * velocity_cfl) / (max(shape)-1)
+                        
+                        napplies = 1
+                        simu = Simulation(start=0.0, end=napplies*dt, dt0=dt)
+                        simu.initialize()
+                        while not simu.is_over:
+                            split.apply(simulation=simu)
+                            simu.advance()
+                        
+                        offsets = (-Vi*napplies*dt).tolist()
+                        dsref.initialize(self.__scalar_init, offsets=offsets)
+
+                        #print 'FINAL', dsout[0][view].sum()
+                        #print 'EXPECTED', dsref[0][view].sum()
+                        print Vi, dsout.distance(dsref)
+                        #print
+                        
+                        simu.finalize()
+                        split.finalize()
+                    except:
+                        sys.stdout.write('\bx\n\n')
+                        sys.stdout.flush()
+                        raise
     
     @classmethod
     def _check_output(cls, impl, op, refin_buffers, refout_buffers, out_buffers):
@@ -189,30 +215,34 @@ class TestDirectionalAdvectionOperator(object):
         self._test(dim=1, is_inplace=False)
     def test_advec_2D_out_of_place(self):
         self._test(dim=2, is_inplace=False)
-    def test_advec_2D_out_of_place(self):
+    def test_advec_3D_out_of_place(self):
         self._test(dim=3, is_inplace=False)
     
     def test_advec_1D_inplace(self):
         self._test(dim=1, is_inplace=True)
     def test_advec_2D_inplace(self):
         self._test(dim=2, is_inplace=True)
-    def test_advec_2D_inplace(self):
+    def test_advec_3D_inplace(self):
         self._test(dim=3, is_inplace=True)
     
     def perform_tests(self):
-        self.test_advec_1D_out_of_place()
-        self.test_advec_2D_out_of_place()
-        self.test_advec_3D_out_of_place()
+        #self.test_advec_1D_out_of_place()
+        #self.test_advec_2D_out_of_place()
+        #self.test_advec_3D_out_of_place()
         
-        self.test_advec_1D_inplace()
+        #self.test_advec_1D_inplace()
         self.test_advec_2D_inplace()
-        self.test_advec_3D_inplace()
+        #self.test_advec_3D_inplace()
     
 if __name__ == '__main__':
     TestDirectionalAdvectionOperator.setup_class(enable_extra_tests=False, 
                                       enable_debug_mode=False)
     
     test = TestDirectionalAdvectionOperator()
-    test.perform_tests()
+    
+    with printoptions(threshold=10000, linewidth=240, 
+                      nanstr='nan', infstr='inf', 
+                      formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
+        test.perform_tests()
 
     TestDirectionalAdvectionOperator.teardown_class()
diff --git a/hysop/parameters/numpy_parameter.py b/hysop/parameters/numpy_parameter.py
index b0d73aaf0..045a47e15 100644
--- a/hysop/parameters/numpy_parameter.py
+++ b/hysop/parameters/numpy_parameter.py
@@ -93,17 +93,17 @@ class NumpyParameter(Parameter):
                 parameter_types=parameter_types, allow_None=False,
                 initial_value=initial_value, **kwds)
 
-        if obj.obj_initialized:
+        #if obj.obj_initialized:
             # object was already created, check shape, dtype, min and max values
-            assert obj.shape == shape, 'shape mismatch.'
-            assert obj.dtype == dtype, 'dtype mismatch.'
-            assert obj.min_value == min_value, 'min_value mismatch.'
-            assert obj.max_value == max_value, 'max_value mismatch.'
-            assert obj.ignore_nans == ignore_nans, 'ignore_nans mismatch.'
-        else:
-            obj._min_value = min_value
-            obj._max_value = max_value
-            obj._ignore_nans = ignore_nans
+            #assert obj.shape == shape, 'shape mismatch.'
+            #assert obj.dtype == dtype, 'dtype mismatch.'
+            #assert obj.min_value == min_value, 'min_value mismatch.'
+            #assert obj.max_value == max_value, 'max_value mismatch.'
+            #assert obj.ignore_nans == ignore_nans, 'ignore_nans mismatch.'
+        #else:
+        obj._min_value = min_value
+        obj._max_value = max_value
+        obj._ignore_nans = ignore_nans
 
         return obj
 
diff --git a/hysop/parameters/parameter.py b/hysop/parameters/parameter.py
index bd0005053..08d07134c 100644
--- a/hysop/parameters/parameter.py
+++ b/hysop/parameters/parameter.py
@@ -6,10 +6,10 @@ Parameters description.
 from abc import ABCMeta, abstractmethod
 from hysop import dprint, vprint
 from hysop.tools.types import check_instance, to_tuple
-from hysop.tools.handle import RegisteredObject
+from hysop.tools.handle import TaggedObject
 from hysop.tools.variable import Variable, VariableTag
 
-class Parameter(RegisteredObject, VariableTag):
+class Parameter(TaggedObject, VariableTag):
     """
     A parameter is a value of a given type that may change value as simulation advances.
     Parameters are only available on the host backend.
@@ -18,7 +18,7 @@ class Parameter(RegisteredObject, VariableTag):
 
     def __new__(cls, name, parameter_types, 
                     initial_value=None, allow_None=False,
-                    ensure_unique=True, **kwds):
+                    **kwds):
         """
         Create or get an existing Parameter with a specific name
         and type.
@@ -63,28 +63,35 @@ class Parameter(RegisteredObject, VariableTag):
 
         parameter_types = tuple(set(parameter_types))
 
-        obj = super(Parameter, cls).__new__(cls, name=name, 
+        obj = super(Parameter, cls).__new__(cls, 
                 tag_prefix='p', variable_kind=Variable.PARAMETER, **kwds)
         
-        if obj.obj_initialized:
+        #if obj.obj_initialized:
             # object was already created, check name and initial value
             # unless ensure_unique has been set.
-            assert (obj.name == name), 'FATAL ERROR: registered object name mismatch.'
-            if ensure_unique:
-                msg='A Parameter with the same name {} has already been created, '
-                msg+='and ensure_unique flag has been set.'
-                msg=msg.format(name)
-                raise RuntimeError(msg)
-            else:
-                assert (obj._parameter_types == parameter_types), 'Allowed types mismatch.'
-                assert (initial_value is None) or (obj._value == initial_value), 'Initial value does not match old one.'
-        else:
-            obj._name = name
-            obj._parameter_types = parameter_types
-            obj._value = initial_value
+            #assert (obj.name == name), 'FATAL ERROR: registered object name mismatch.'
+            #if ensure_unique:
+                #msg='A Parameter with the same name {} has already been created, '
+                #msg+='and ensure_unique flag has been set.'
+                #msg=msg.format(name)
+                #raise RuntimeError(msg)
+            #else:
+                #assert (obj._parameter_types == parameter_types), 'Allowed types mismatch.'
+                #assert (initial_value is None) or (obj._value == initial_value), 'Initial value does not match old one.'
+        #else:
+        obj._name = name
+        obj._parameter_types = parameter_types
+        obj._value = initial_value
 
         return obj
 
+    def __eq__(self, other):
+        return (self is other)
+    def __neq__(self, other):
+        return (self is not other)
+    def __hash__(self):
+        return id(self)
+
     def set_value(self, value):
         """Set the value of this Parameter object."""
         if not isinstance(value, self._parameter_types):
diff --git a/hysop/simulation.py b/hysop/simulation.py
index e6a2e8971..0f1421471 100644
--- a/hysop/simulation.py
+++ b/hysop/simulation.py
@@ -33,8 +33,9 @@ from hysop import dprint
 from hysop.deps import sys
 from hysop.constants import HYSOP_REAL
 from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.tools.numpywrappers import npw
 
-eps = sys.float_info.epsilon
+eps = npw.finfo(HYSOP_REAL).eps
 """Machine epsilon, used to compare times."""
 
 
@@ -115,7 +116,7 @@ class Simulation(object):
         self._dt0 = dt0
 
         dt_name = '{}_dt'.format(name) if (name is not None) else 'dt'
-        self.dt = ScalarParameter(name=dt_name, dtype=HYSOP_REAL, min_value=1e-12,
+        self.dt = ScalarParameter(name=dt_name, dtype=HYSOP_REAL, min_value=eps,
                 initial_value=dt0, ensure_unique=kwds.pop('ensure_unique', True))
         self.name = name
         
@@ -199,11 +200,16 @@ class Simulation(object):
         self.update_time_step(self._dt0)
         self.tkp1 = self.start + self.time_step
         assert self.tkp1 <= self.end
+        if abs(self._dt0 - self.end) <= self.tol:
+            self._next_is_last = True
+        else:
+            self._next_is_last = False
+        
         self.time = self.tkp1
         self.is_over = False
         self.current_iteration = 0
         self._is_ready = True
-        self._next_is_last = False
+
 
     def finalize(self):
         """Use this function when you need to call an hdf i/o operator
diff --git a/hysop/tools/transposition_states.py b/hysop/tools/transposition_states.py
index ed3ab0df8..153ae6dc6 100644
--- a/hysop/tools/transposition_states.py
+++ b/hysop/tools/transposition_states.py
@@ -143,10 +143,9 @@ class TranspositionState(object):
         return hash(self._axes)
     
     def __str__(self):
-        return 'TranspositionState{}D.{}'.format(self.dim, 
-                ''.join(DirectionLabels[i] for i in self.axes))
+        return ''.join(DirectionLabels[i] for i in self.axes)
     def __repr__(self):
-        return self.__str__()
+        return 'TranspositionState{}D.{}'.format(self.dim, self.__str__())
 
     axes = property(_get_axes)
     dim  = property(_get_dim)
diff --git a/hysop/topology/cartesian_topology.py b/hysop/topology/cartesian_topology.py
index 026ed5034..450ba1c46 100644
--- a/hysop/topology/cartesian_topology.py
+++ b/hysop/topology/cartesian_topology.py
@@ -171,7 +171,9 @@ class CartesianTopologyState(TopologyState):
     axes   = property(_get_axes, _set_axes)
     basis  = property(_get_basis, _set_basis)
     dim    = property(_get_dim)
+
     tstate = property(_get_tstate)
+    transposition_state = property(_get_tstate)
 
 class CartesianTopologyView(TopologyView):
     """
-- 
GitLab