From cfc5f1cfe2b8e37f24747a8ab27f21ba3dfd3aa9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Mon, 13 May 2013 14:07:40 +0000
Subject: [PATCH] Update analytical fields and operators
---
HySoP/hysop/fields/__init__.py | 22 ++++++++-
HySoP/hysop/fields/analytical.py | 56 -----------------------
HySoP/hysop/fields/continuous.py | 43 +++++++++++++++--
HySoP/hysop/fields/discrete.py | 2 +-
HySoP/hysop/fields/scalar.py | 18 ++++----
HySoP/hysop/fields/vector.py | 18 ++++----
HySoP/hysop/operator/analytic.py | 23 ++++++++--
HySoP/hysop/operator/discrete/analytic.py | 13 ++----
8 files changed, 98 insertions(+), 97 deletions(-)
delete mode 100644 HySoP/hysop/fields/analytical.py
diff --git a/HySoP/hysop/fields/__init__.py b/HySoP/hysop/fields/__init__.py
index 16ef1de56..d709fa847 100644
--- a/HySoP/hysop/fields/__init__.py
+++ b/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
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
deleted file mode 100644
index 98a09a4b9..000000000
--- a/HySoP/hysop/fields/analytical.py
+++ /dev/null
@@ -1,56 +0,0 @@
-"""
-@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__
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index d51ef693f..1da1a1f79 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/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:
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
index a7884643d..1ff8f84cb 100644
--- a/HySoP/hysop/fields/discrete.py
+++ b/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):
diff --git a/HySoP/hysop/fields/scalar.py b/HySoP/hysop/fields/scalar.py
index 0383f4755..4be9aaf69 100644
--- a/HySoP/hysop/fields/scalar.py
+++ b/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):
"""
diff --git a/HySoP/hysop/fields/vector.py b/HySoP/hysop/fields/vector.py
index 8c4fd90f6..e4b4c7210 100644
--- a/HySoP/hysop/fields/vector.py
+++ b/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):
"""
diff --git a/HySoP/hysop/operator/analytic.py b/HySoP/hysop/operator/analytic.py
index d5a0b3ccf..e7c23ea0b 100644
--- a/HySoP/hysop/operator/analytic.py
+++ b/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,
diff --git a/HySoP/hysop/operator/discrete/analytic.py b/HySoP/hysop/operator/discrete/analytic.py
index d663652b8..4d729d98c 100644
--- a/HySoP/hysop/operator/discrete/analytic.py
+++ b/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:
--
GitLab