external force

examples/sediment_deposit/sediment_deposit.py

@@ -248,7 +248,7 @@ def compute(args):
     g = 1.0
     Fext = SymbolicExternalForce(name='S', Fext=(0,-g*Ss),
                                    diffusion = {S: nu_S})
-    external_force = SpectralExternalForce(vorticity=Ws, dt=dt,
+    external_force = SpectralExternalForce(vorticity=vorti, dt=dt,
                                    Fext=Fext, Finf=True,
                                    implementation=impl,
                                    variables={vorti: npts, S: npts},

hysop/backend/device/opencl/operator/external_force.py

@@ -4,7 +4,7 @@ from hysop.tools.numerics import float_to_complex_dtype
 from hysop.tools.numpywrappers import npw
 from hysop.tools.decorators  import debug
 from hysop.tools.warning import HysopWarning
-from hysop.operator.base.curl import SpectralExternalForceOperatorBase
+from hysop.operator.base.external_force import SpectralExternalForceOperatorBase
 from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
 from hysop.core.graph.graph import op_apply
 from hysop.core.memory.memory_request import MemoryRequest
@@ -22,82 +22,13 @@ class OpenClSpectralExternalForce(SpectralExternalForceOperatorBase, OpenClSymbo
     
     @debug
     def __init__(self, **kwds):
-        """
-        What is computed in 2D:
-            FEXT_X = fft(Fext_X)
-            tmp    = dFEXT/dy
-            FEXT_Y = fft(Fext_Y)
-            tmp   -= dFEXT/dx
-            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
-            Wz  += dt*tmp
-        What is computed in 3D:
-            FEXT_X = fft(Fext_X)
-            FEXT_Y = fft(Fext_Y)
-            FEXT_Z = fft(Fext_Z)
-            
-            tmp = 
-            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
-            Wx += dt*tmp
-
-            tmp = 
-            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
-            Wy += dt*tmp
-
-            tmp = 
-            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
-            Wz += dt*tmp
-        """
         super(OpenClSpectralExternalForce, self).__init__(**kwds)
 
-        assert (len(self.forward_transforms) % 2 == 0)
-        N = len(self.forward_transforms)//2
-        assert len(self.K)==2*N
-
-        kernel_names = ()
-        for (i,(Ft,(tg,Ki))) in enumerate(zip(self.forward_transforms, self.K)):
-            Fhs = Ft.output_symbolic_array('F{}_hat'.format(i))
-            
-            kname = 'filter_curl_{}d_{}'.format(Fhs.dim, i)
-            kernel_names += (kname,)
-            
-            is_complex = Ki.is_complex
-            Ki = tg._indexed_wave_numbers[Ki]
-
-            if is_complex:
-                expr = ComplexMul(Ki, Fhs)
-            else:
-                expr = Ki*Fhs
-
-            if (i<N):
-                expr = Assignment(Fhs, +expr)
-            else:
-                expr = Assignment(Fhs, -expr)
-
-            self.require_symbolic_kernel(kname, expr)
-
-        self._kernel_names = kernel_names
-    
     @debug
     def setup(self, work):
         super(OpenClSpectralExternalForce, self).setup(work)
-        curl_filters = ()
-        for kname in self._kernel_names:
-            kernel, _ = self.symbolic_kernels[kname]
-            kernel = functools.partial(kernel, queue=self.cl_env.default_queue)
-            curl_filters += (kernel,)
-        self.curl_filters = curl_filters
-        self.exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
-        assert len(self.forward_transforms)==len(self.backward_transforms)==len(curl_filters)
 
     @op_apply
     def apply(self, **kwds):
-        """Solve the ExternalForce equation."""
         super(OpenClSpectralExternalForce, self).apply(**kwds)
-        for (Ft,filter_curl,Bt) in zip(self.forward_transforms,
-                                       self.curl_filters,
-                                       self.backward_transforms):
-            evt = Ft()
-            evt = filter_curl()
-            evt = Bt()
-        if (self.exchange_ghosts is not None):
-            evt = self.exchange_ghosts()
+

hysop/operator/adapt_timestep.py

@@ -76,8 +76,7 @@ class TimestepCriteria(ComputationalGraphOperator):
 class ConstantTimestepCriteria(TimestepCriteria):
 
     @debug
-    def __init__(self, cst, parameter, criteria,
-                    Finf=None,
+    def __init__(self, cst, parameter, Finf,
                     name=None, pretty_name=None, **kwds):
         """
         Initialize a ConstantTimestepCriteria.
@@ -98,8 +97,10 @@ class ConstantTimestepCriteria(TimestepCriteria):
         kwds: dict
             Base class arguments.
         """
+        if isinstance(cst, (int, long)):
+            cst = float(cst)
         assert (cst > 0.0), 'negative cst factor.'
-        check_instance(cst, (float, np.ndarray, list, tuple))
+        check_instance(cst, (float, npw.ndarray, list, tuple))
         check_instance(Finf, (ScalarParameter, TensorParameter))
         if isinstance(Finf, ScalarParameter):
             assert isinstance(cst, float)
@@ -107,19 +108,18 @@ class ConstantTimestepCriteria(TimestepCriteria):
         else:
             is_scalar = False
             if isinstance(cst, float):
-                cst = np.full(shape=Finf.shape, dtype=Finf.dtype, fill_value=cst)
+                cst = npw.full(shape=Finf.shape, dtype=Finf.dtype, fill_value=cst)
             if isinstance(cst, (list, tuple)):
                 assert Finf.ndim == 1
-                cst = np.asarray(cst)
+                cst = npw.asarray(cst)
             msg='Shape mismatch between parameter {} and cst {}.'
             msg=msg.format(Finf.shape, cst.shape)
             assert Finf.shape == cst.shape, msg
 
-        input_params[Finf.name] = Finf
         name = first_not_None(name, 'CST')
         pretty_name = first_not_None(pretty_name, name)
         super(ConstantTimestepCriteria, self).__init__(name=name, pretty_name=pretty_name,
-                                                 input_params=input_params,
+                                                 input_params={Finf.name: Finf},
                                                  output_params={parameter.name: parameter},
                                                  parameter=parameter, **kwds)
         self.cst = cst
@@ -133,13 +133,13 @@ class ConstantTimestepCriteria(TimestepCriteria):
         if self.is_scalar:
             assert Finf >= 0
             if (Finf == 0):
-                return np.inf
+                return npw.inf
             else:
                 return cst/Finf
         else:
             assert Finf.min() >= 0
             mask = (Finf!=0)
-            dt = np.zeros_like(cst, fill_value=np.inf)
+            dt = npw.zeros_like(cst, fill_value=npw.inf)
             dt[mask] = cst[mask] / Finf[mask]
             return dt.min()
 
@@ -511,8 +511,7 @@ class AdaptiveTimeStep(ComputationalGraphNodeGenerator):
                 name=param_name, pretty_name=param_pretty_name, 
                 basename=name)
         criteria = ConstantTimestepCriteria(cst=cst, Finf=Finf,
-            parameter=parameter,
-            name=name, pretty_name=pretty_name, **kwds)
+            parameter=parameter, name=name, pretty_name=pretty_name, **kwds)
         self._push_criteria(parameter.name, criteria)
 
     def push_cfl_criteria(self, cfl, Fmin=None, Fmax=None, Finf=None,

hysop/operator/base/min_max.py

@@ -27,7 +27,7 @@ class MinMaxFieldStatisticsBase(object):
 
     @classmethod
     def build_parameters(cls, field, components, all_quiet,
-                    Fmin, Fmax, Finf, pbasename, ppbasename):
+                    Fmin, Fmax, Finf, pbasename, ppbasename, dtype=None):
         if (    ((Fmin is None) or (Fmin is False)) 
             and ((Fmax is None) or (Fmax is False))
             and ((Finf is None) or (Finf is False))):
@@ -37,8 +37,10 @@ class MinMaxFieldStatisticsBase(object):
             msg+=' tensor parameter.'
             raise ValueError(msg)
         
-        pbasename  = first_not_None(pbasename,  field.name)
-        ppbasename = first_not_None(ppbasename, field.pretty_name.decode('utf-8'))
+        if (field is not None):
+            dtype      = first_not_None(dtype,      field.dtype)
+            pbasename  = first_not_None(pbasename,  field.name)
+            ppbasename = first_not_None(ppbasename, field.pretty_name.decode('utf-8'))
         
         def make_param(k, quiet):
             return TensorParameter(name=names[k], pretty_name=pretty_names[k],
@@ -56,7 +58,9 @@ class MinMaxFieldStatisticsBase(object):
             'Finf': u'|{}|\u208a'.format(ppbasename),
         }
         
-        components = to_tuple(first_not_None(components, range(field.nb_components)))
+        if (field is not None):
+            components = first_not_None(components, range(field.nb_components))
+        components = to_tuple(components)
         nb_components = len(components)
 
         parameters = {}

hysop/operator/external_force.py

 
+import sympy as sm
 from abc import ABCMeta, abstractmethod
-from hysop.fields
-from hysop.symbolic.field import AppliedSymbolicField, SymbolicDiscreteField
+from hysop.constants import Implementation
+from hysop.fields.continuous_field import Field
+from hysop.symbolic.field import AppliedSymbolicField, SymbolicField
 from hysop.symbolic.parameter import SymbolicScalarParameter
+from hysop.parameters.tensor_parameter import TensorParameter
 from hysop.parameters.scalar_parameter import ScalarParameter
 from hysop.operator.min_max import MinMaxFieldStatisticsBase
+from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.tools.types import first_not_None, to_tuple, check_instance
+from hysop.tools.sympy_utils import nabla
+from hysop.tools.interface import NamedObjectI
+from hysop.tools.decorators import debug
 
 class ExternalForce(NamedObjectI):
     __metaclass__ = ABCMeta
@@ -16,18 +25,12 @@ class ExternalForce(NamedObjectI):
         self._dim  = dim
         self._Fext = Fext
 
-        for f in self.input_fields:
+        for f in self.input_fields():
             msg='Expected field dimension to be {} but got {} from field {}.'
-            msg=msg.format(dim, field.dim, field.short_description())
+            msg=msg.format(dim, f.dim, f.short_description())
             assert f.dim == dim, msg
 
     
-    @property
-    def name(self):
-        return self._diffusion
-    @property
-    def pretty_name(self):
-        return self._diffusion
     @property
     def Fext(self):
         return self._Fext
@@ -52,18 +55,18 @@ class ExternalForce(NamedObjectI):
         pass
 
 
-class SymbolicExternalForce(object):
+class SymbolicExternalForce(ExternalForce):
     def __init__(self, name, Fext, diffusion=None, **kwds):
         """
         Specify an external force as a tuple of symbolic expressions.
 
         2D ExternalForce example:  
             1) Fext = -rho*g*e_y where rho is a field and g a constant
-                Fext = (0, -rho.s*g)
+                Fext = (0, -rho.s()*g)
             2) Fext = (Rs*S+C)*e_y
-                Fext = (0, -Rs*S.s+C.s)
+                Fext = (0, -Rs*S.s()+C.s())
         """
-        check_instance(Fext, tuple)
+        Fext = to_tuple(Fext)
         dim = len(Fext)
 
         Fext = list(Fext)
@@ -72,14 +75,15 @@ class SymbolicExternalForce(object):
                 Fext[i] = 0 # curl(const) = 0
             msg='Expression "{}" contains a SymbolicField, did you forget to apply it ?'
             msg=msg.format(e)
-            assert not e.atoms(SymbolicField), msg
+            if isinstance(e, sm.Basic):
+                assert not e.atoms(SymbolicField), msg
         Fext = tuple(Fext)
 
         super(SymbolicExternalForce, self).__init__(name=name, dim=dim, Fext=Fext, **kwds)
         
         diffusion = first_not_None(diffusion, {})
         for (k,v) in diffusion.iteritems():
-            assert k in self.input_fields
+            assert k in self.input_fields(), k.short_description()
             assert isinstance(v, ScalarParameter)
         self._diffusion = diffusion
 
@@ -95,11 +99,44 @@ class SymbolicExternalForce(object):
         return set()
 
     def _extract_objects(self, obj_type):
-        objs = ()
+        objs = set()
         for e in self.Fext:
-            objs += e.atoms(obj_type)
+            try:
+                objs.update(e.atoms(obj_type))
+            except AttributeError:
+                pass
         return objs
 
+    def short_description(self):
+        return 'SymbolicExternalForce[name={}]'.format(self.name)
+    
+    def long_description(self):
+        sep = '\n      *'
+        expressions = sep + sep.join('F{} = {}'.format(x,e) for (x,e) in zip('xyz',self.Fext))
+        diffusion = sep + sep.join('{}: {}'.format(f.pretty_name, p.pretty_name) 
+                                        for (f,p) in self.diffusion.iteritems())
+        input_fields  = ', '.join(f.pretty_name for f in self.input_fields())
+        output_fields = ', '.join(f.pretty_name for f in self.output_fields())
+        input_params  = ', '.join(p.pretty_name for p in self.input_params())
+        output_params = ', '.join(p.pretty_name for p in self.output_params())
+        
+        ss = \
+        '''SymbolicExternalForce:
+    name:          {}
+    pretty_name:   {}
+    expressions:   {}
+    diffusion:     {}
+    -----------------
+    input_fields:  {}
+    output_fields: {}
+    input_params:  {}
+    output_params: {}
+        '''.format(self.name, self.pretty_name, 
+                expressions, diffusion, 
+                input_fields, output_fields,
+                input_params, output_params)
+        return ss
+
 
 class SpectralExternalForce(ComputationalGraphNodeFrontend):
     """
@@ -120,10 +157,11 @@ class SpectralExternalForce(ComputationalGraphNodeFrontend):
 
     @debug
     def __init__(self, vorticity, Fext, dt, variables,
-                    Fmin=None, Fmax=None, Finf=None,
+                    Fmin=None, Fmax=None, Finf=None, 
+                    all_quiet=False, pbasename=None, ppbasename=None,
                     implementation=None, base_kwds=None, **kwds):
         """
-        Create an operator that computes the curl of a given input field Fext or a set of symbolic expressions.
+        Create an operator that computes the curl of a given input force field Fext.
 
         Only the following configurations are supported:
                      dim   nb_components  |   dim   nb_components
@@ -145,10 +183,15 @@ class SpectralExternalForce(ComputationalGraphNodeFrontend):
             All contained field have to live on the same domain.
         Fext: hysop.operator.external_force.ExternalForce
             Expression of the external force.
-        F...: TensorParameter or boolean, optional
+        dt: ScalarParameter
+            Timestep paramater.
+        F...: ScalarParameter, TensorParameter or boolean, optional
             TensorParameters should match the shape of tmp (see Notes).
             If set to True, the TensorParameter will be generated automatically.
-            Autogenerated TensorParameters that are not required by the user are set to be quiet.
+        all_quiet: bool
+            Force all autogenerated TensorParameter to be quiet.
+            By default, only the autogenerated TensorParameters that are not required 
+            by the user are set to be quiet.
         pbasename: str, optional
             Parameters basename for created parameters.
             Defaults to 'curl_{}'.format(Fext.name).
@@ -165,28 +208,42 @@ class SpectralExternalForce(ComputationalGraphNodeFrontend):
 
         Notes
         -----
-        If Fext.dim == 2, Fext
+        If dim == 2, it is expected that:
+            vorticity has only one component
+            Fext has 2 components
+            Expected parameters are ScalarParameters or TensorParameters of shape (1,)
+        Else if dim == 3:
+            vorticity has 3 components
+            Fext has 3 components
+            Expected parameters are TensorParameter of shape (3,)
         """
         base_kwds = first_not_None(base_kwds, {})
         check_instance(vorticity, Field)
         check_instance(Fext, ExternalForce)
+        check_instance(dt, ScalarParameter)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         
         # Pregenerate parameters so that we can directly store them in self.
         default_pbasename  = 'curl_{}'.format(Fext.name)
-        default_ppbasename = '{}_{}'.format(nabla, Fext.pretty_name)
+        default_ppbasename = u'{}x{}'.format(nabla, Fext.pretty_name)
         pbasename  = first_not_None(pbasename,  default_pbasename)
         ppbasename = first_not_None(ppbasename, default_ppbasename)
-        parameters = MinMaxFieldStatisticsBase.build_parameters(field=field, 
-                components=components, all_quiet=all_quiet,
+        parameters = MinMaxFieldStatisticsBase.build_parameters(field=vorticity, 
+                components=None, dtype=None, all_quiet=all_quiet, 
                 Fmin=Fmin, Fmax=Fmax, Finf=Finf, 
                 pbasename=pbasename, ppbasename=ppbasename)
 
         (Fmin, Fmax, Finf) = tuple(parameters[k] for k in ('Fmin', 'Fmax', 'Finf'))
-
-        super(Curl, self).__init__(Fin=Fin, Fout=Fout, variables=variables,
-                candidate_input_tensors=(Fin,),
-                candidate_output_tensors=(Fout,),
+        
+        check_instance(Fmin, TensorParameter, allow_none=True)
+        check_instance(Fmax, TensorParameter, allow_none=True)
+        check_instance(Finf, TensorParameter, allow_none=True)
+
+        super(SpectralExternalForce, self).__init__(vorticity=vorticity, 
+                Fext=Fext, dt=dt, variables=variables,
+                Fmin=Fmin, Fmax=Fmax, Finf=Finf,
+                candidate_input_tensors=(vorticity,),
+                candidate_output_tensors=(vorticity,),
                 implementation=implementation, 
                 base_kwds=base_kwds, **kwds)
 

hysop/tools/interface.py

@@ -39,17 +39,18 @@ class NamedObjectI(object):
                     latex_name=latex_name, var_name=var_name)
         return obj
     
-    def rename(self, name, pretty_name=None, latex_name=None):
+    def rename(self, name, pretty_name=None, latex_name=None, var_name=None):
         """Change the names of this object."""
         check_instance(name, str)
         check_instance(pretty_name, (str,unicode), allow_none=True)
         check_instance(latex_name, str, allow_none=True)
-        if isinstance(pretty_name, unicode):
-            pretty_name = pretty_name.encode('utf-8')
-        check_instance(pretty_name, str)
         
         pretty_name = first_not_None(pretty_name, name)
         latex_name  = first_not_None(latex_name, name)
+        
+        if isinstance(pretty_name, unicode):
+            pretty_name = pretty_name.encode('utf-8')
+        check_instance(pretty_name, str)
 
         self._name   = name
         self._pretty_name = pretty_name