Update analytical fields and operators

HySoP/hysop/fields/__init__.py

@@ -19,7 +19,6 @@
 # \code
 # from parmepy.domain.box import Box
 # from parmepy.fields.continuous import Field
-# from parmepy.analytical import AnalyticalField
 #
 # # First, the domain : a 3D box
 # dom = Box()
@@ -35,11 +34,30 @@
 #     valz = sin(z)
 #     return valx, valy, valz
 #
-# vorticity = AnalyticalField(dom, isVector=True, name='reference,
+# vorticity = Field(dom, isVector=True, name='reference,
 #     formula=compute_ref)
 #
 # # Compute the value of vorticity at point of coordinates (1,2,0) :
 # wx, wy, wz = vorticity.value(1.,2.,0.)
+#
+# # Modify the formula and add some extra arguments:
+# 
+# def compute_ref2(x,y,z,t,dt):
+#     valx = sin(x)*t
+#     valy = sin(y) + dt
+#     valz = sin(z)
+#     return valx, valy, valz
+#
+# vorticity.setFormula(compute_ref2)
+#
+# t = 4
+# dt = 8.3
+# # Compute vorticity using compute_ref2 for all the available
+# # discretizations.
+# # Warning : space coordinates are implicit but all extra-parameters
+# # must be explicitly given, in the same order as defined in compute_ref2.
+# # No keywords arguments allowed.
+# vorticity.initialize(t,dt)
 # \endcode
 #
 # Up to now, these are just descriptions of the fields. No discretisations, no

HySoP/hysop/fields/analytical.py

-"""
-@file fields/analytical.py
-
-Continuous variable description defined with an analytic formula.
-"""
-from parmepy.constants import debug
-from parmepy.fields.continuous import Field
-
-
-class AnalyticalField(Field):
-    """
-    A continuous variable defined with an analytic formula.
-
-    """
-
-    @debug
-    def __init__(self, domain, formula, name=None, isVector=False):
-        """
-        Create an AnalyticalField from a formula.
-
-        @param domain : variable application domain.
-        @param formula : a user-defined python function
-        @param name : name of the variable (may be used for vtk output).
-        @param isVector : true if the field is a vector.
-
-        @note formula is used in ScalarField or VectorField as a vectorized
-        function by numpy.
-        """
-        Field.__init__(self, domain, name, isVector)
-        ## Analytic formula.
-        self.formula = formula
-
-    def value(self, *pos):
-        """
-        Evaluation of the variable at a given position.
-
-        @param pos : Position to evaluate.
-        @return formula(pos): value of the formula at the given position.
-        """
-        return self.formula(*pos)
-
-    @debug
-    def initialize(self):
-        for dF in self.discreteFields.values():
-            dF.initialize(self.formula)
-
-    def __str__(self):
-        """ToString method"""
-        s = " --- Analytical variable --- \n"
-        s += Field.__str__(self)
-        return s
-
-if __name__ == "__main__":
-    print __doc__
-    print "- Provided class : AnalyticalField"
-    print AnalyticalField.__doc__

HySoP/hysop/fields/continuous.py

@@ -38,12 +38,14 @@ class Field(object):
         return object.__new__(cls, *args, **kw)
 
     @debug
-    def __init__(self, domain, name=None, isVector=False):
+    def __init__(self, domain, formula=None, name=None, isVector=False):
         """
         Create a scalar of vector field which dimension is determined
         by its domain of definition.
 
         @param domain : physical domain where this field is defined.
+        @param formula : a user-defined python function
+        (default = 0 all over the domain)
         @param name : name of the variable.
         @param isVector : true if the field is a vector.
         """
@@ -63,6 +65,8 @@ class Field(object):
         self.discreteFields = {}
         ## Is this field a vector field?
         self.isVector = isVector
+        ## Analytic formula.
+        self._formula = formula
 
     @debug
     def discretize(self, topo):
@@ -88,22 +92,51 @@ class Field(object):
 
         return self.discreteFields[topo]
 
+    def setFormula(self, formula):
+        """
+        Set (or reset) a user-defined function to compute field
+        values.
+        @param formula : a user-defined python function
+        Warning: if set, this formula will
+        overwrite the (optional) formula given during
+        construction of this field.
+
+        """
+        self._formula = formula
+
     @debug
-    def initialize(self):
+    def initialize(self, *args):
         """
         Initialize all the discrete fields associated to this continuous field.
+        using the formula set during construction or with setFormula method.
+        If formula is not set, field values are set to zero.
 
-        @see VectorField.initialize
-        @see ScalarField.initialize
+        @param *args : see parmepy.fields
         """
         for dF in self.discreteFields.values():
-            dF.initialize()
+            dF.initialize(self._formula, *args)
+
+    def value(self, *pos):
+        """
+        Evaluation of the field at a given position, according
+        to the analytic formula given during construction.
+
+        @param pos : Position to evaluate.
+        @return formula(pos): value of the formula at the given position.
+        """
+        if self._formula is not None:
+            return self._formula(*pos)
+        else:
+            return 0.0
 
     def __str__(self):
         """ Field info display """
 
         s = '[' + str(main_rank) + '] " ' + self.name + ' ", '
         s += str(self.dimension) + 'D '
+        if self._formula is not None:
+            s += 'an analytic formula is associated to this field:'
+            s += str(self._formula.func_name)
         if self.isVector:
             s += 'vector field '
         else:

HySoP/hysop/fields/discrete.py

@@ -103,7 +103,7 @@ class DiscreteField(object):
 
     @debug
     @abstractmethod
-    def initialize(self, formula=None):
+    def initialize(self, formula=None, *args, **kwargs):
         """Initialize values with given formula."""
 
     def get_data_method(self):

HySoP/hysop/fields/scalar.py

@@ -46,7 +46,7 @@ class ScalarField(DiscreteField):
         self.data.__setitem__(i, value)
 
     @debug
-    def initialize(self, formula=None):
+    def initialize(self, formula=None, *args):
         """
         Initialize values with given formula.
 
@@ -58,22 +58,20 @@ class ScalarField(DiscreteField):
         @li 2D : float formula(float x, float y)
         where x,y,z are point coordinates.
         """
-        if formula is not None and not self.contains_data:
-            if __VERBOSE__:
-                print "...",
+        if formula is not None:
             v_formula = np.vectorize(formula)
             if self.dimension == 3:
-                self.data[...] = v_formula(*self.topology.mesh.coords)
+                self.data[...] = v_formula(*(self.topology.mesh.coords + args))
             elif self.dimension == 2:
-                self.data[...] = v_formula(*self.topology.mesh.coords)
+                self.data[...] = v_formula(*(self.topology.mesh.coords + args))
             else:
-                self.data[...] = formula(self.topology.mesh.coords)
+                self.data[...] = formula(*(self.topology.mesh.coords + args))
             self.data[...] = np.asarray(self.data,
                                         dtype=PARMES_REAL, order=ORDER)
-            self.contains_data = True
         else:
-            if __VERBOSE__:
-                print "No formula",
+            self.data[...] = 0.0
+
+        self.contains_data = True
 
     def norm(self):
         """

HySoP/hysop/fields/vector.py

@@ -3,7 +3,6 @@
 
 Discrete vector field.
 """
-from parmepy import __VERBOSE__
 from parmepy.constants import np, ORDER, PARMES_REAL, debug
 from parmepy.fields.discrete import DiscreteField
 from numpy import linalg as la
@@ -77,7 +76,7 @@ class VectorField(DiscreteField):
             self.data[i[0]].__setitem__(tuple(i[1:]), value)
 
     @debug
-    def initialize(self, formula=None):
+    def initialize(self, formula=None, *args):
         """
         Initialize values with given formula.
 
@@ -89,13 +88,11 @@ class VectorField(DiscreteField):
         @li 2D : float, float formula(float x, float y)
         where x,y,z are point coordinates.
         """
-        if formula is not None and not self.contains_data:
-            if __VERBOSE__:
-                print "...",
+        if formula is not None:
             v_formula = np.vectorize(formula)
             if self.dimension == 3:
                 self.data[0][...], self.data[1][...], self.data[2][...] = \
-                    v_formula(*self.topology.mesh.coords)
+                    v_formula(*(self.topology.mesh.coords + args))
                 self.data[0][...] = np.asarray(self.data[0],
                                                dtype=PARMES_REAL, order=ORDER)
                 self.data[1][...] = np.asarray(self.data[1],
@@ -104,15 +101,16 @@ class VectorField(DiscreteField):
                                                dtype=PARMES_REAL, order=ORDER)
             elif self.dimension == 2:
                 self.data[0][...], self.data[1][...] = \
-                    v_formula(*self.topology.mesh.coords)
+                    v_formula(*(self.topology.mesh.coords + args))
                 self.data[0][...] = np.asarray(self.data[0],
                                                dtype=PARMES_REAL, order=ORDER)
                 self.data[1][...] = np.asarray(self.data[1],
                                                dtype=PARMES_REAL, order=ORDER)
-            self.contains_data = True
         else:
-            if __VERBOSE__:
-                print "No formula",
+            for d in xrange(self.dimension):
+                self.data[d][...] = 0.0
+
+        self.contains_data = True
 
     def norm(self):
         """

HySoP/hysop/operator/analytic.py

@@ -14,22 +14,35 @@ class Analytic(Operator):
     """
 
     @debug
-    def __init__(self, variables, formula, resolutions, method='',
+    def __init__(self, variables, resolutions, formula=None, method=None,
                  **other_config):
         """
         Create an operator using an analytic formula to compute field(s).
         @param formula : the formula to be applied
         @param resolutions : list of resolutions (one per variable)
         """
-        if not callable(formula) and method.find('gpu') == -1:
-            raise ValueError("Analytic formula must be a function")
+
         if isinstance(variables, list):
             Operator.__init__(self, variables, method)
         else:
             Operator.__init__(self, [variables], method)
+
+        if formula is not None:
+            ## A formula applied to all variables of this operator
+            self.formula = formula
+            for v in self.variables:
+                v.setFormula(formula)
+        else:
+            self.formula = self.variables[0]._formula
+            for v in self.variables:
+                assert v._formula is self.formula
+
+        if self.formula is None and self.method.find('gpu') == -1:
+            raise ValueError("A formula must be set.")
+
+        ## Dict of resolutions (one per variable)
         self.resolutions = resolutions
         self.config = other_config
-        self.formula = formula
 
     @debug
     def setUp(self):
@@ -46,7 +59,7 @@ class Analytic(Operator):
             self.discreteFields[v] = v.discretize(topo)
 
         # Select and create discrete operator
-        if self.method.find('gpu') >= 0:
+        if self.method is not None and self.method.find('gpu') >= 0:
             from parmepy.gpu.gpu_analytic import GPUAnalytic
             self.discreteOperator = GPUAnalytic(self.discreteFields.values(),
                                                 method=self.method,

HySoP/hysop/operator/discrete/analytic.py

@@ -31,18 +31,14 @@ class Analytic_D(DiscreteOperator):
         self.vformula = np.vectorize(self.formula)
 
     @debug
-    def apply(self, t=None, dt=None, ite=None):
+    def apply(self, *args):
         cTime = time.time()
         for df in self.variables:
-            args = list(df.topology.mesh.coords)
-            args.append(t)
-            args.append(dt)
-            args.append(ite)
             datatype = df.data[0].dtype
             if df.isVector:
                 if df.topology.dim == 3:
                     df.data[0][...], df.data[1][...], df.data[2][...] = \
-                        self.vformula(*args)
+                        self.vformula(*(df.topology.mesh.coords + args))
                     df.data[0][...] = np.asarray(df.data[0],
                                                  dtype=datatype, order=ORDER)
                     df.data[1][...] = np.asarray(df.data[1],
@@ -50,13 +46,14 @@ class Analytic_D(DiscreteOperator):
                     df.data[2][...] = np.asarray(df.data[2],
                                                  dtype=datatype, order=ORDER)
                 elif df.topology.dim == 2:
-                    df.data[0][...], df.data[1][...] = self.vformula(*args)
+                    df.data[0][...], df.data[1][...] = \
+                        self.vformula(*(df.topology.mesh.coords + args))
                     df.data[0][...] = np.asarray(df.data[0],
                                                  dtype=datatype, order=ORDER)
                     df.data[1][...] = np.asarray(df.data[1],
                                                  dtype=datatype, order=ORDER)
             else:
-                df.data[...] = self.vformula(*args)
+                df.data[...] = self.vformula(*(df.topology.mesh.coords + args))
                 df.data[...] = np.asarray(df.data,
                                           dtype=datatype, order=ORDER)
         for df in self.output: