diff --git a/examples/sediment_deposit/sediment_deposit.py b/examples/sediment_deposit/sediment_deposit.py
index e5f827dc5d294e0ed53dc73bb6fc73ece8e09268..82c66c9a37183e53b6ad8716867bf7bb22148385 100644
--- a/examples/sediment_deposit/sediment_deposit.py
+++ b/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},
diff --git a/hysop/backend/device/opencl/operator/external_force.py b/hysop/backend/device/opencl/operator/external_force.py
index d5e10fb66ecca28ed4351b6693449379bc335b0c..64238517609fb787715242becdee0a8daadaaceb 100644
--- a/hysop/backend/device/opencl/operator/external_force.py
+++ b/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()
+
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 820bd4946bce6c00d8c4afed5f0935bc4e6c5241..68fee1950cedf3c6ed917835882763e1fa596bc3 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/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,
diff --git a/hysop/operator/base/min_max.py b/hysop/operator/base/min_max.py
index 0bbed36412c78c2f184581afe5b517c4810e1e79..964e2fcf83ed6d462ee1e3dca751ce362421cd3a 100644
--- a/hysop/operator/base/min_max.py
+++ b/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 = {}
diff --git a/hysop/operator/external_force.py b/hysop/operator/external_force.py
index 564a113a51001f0e02dd172317afeaadbe120471..bb77ad21c5bfa09be1ba8c136a8ae342c36d7c51 100644
--- a/hysop/operator/external_force.py
+++ b/hysop/operator/external_force.py
@@ -1,10 +1,19 @@
+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)
diff --git a/hysop/tools/interface.py b/hysop/tools/interface.py
index be755838b7aac7819c6e8ceb05f92a4865ea3c34..dcd1334ab0550d5c37fdd01b540812049136b94d 100644
--- a/hysop/tools/interface.py
+++ b/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