diff --git a/examples/sediment_deposit/sediment_deposit.py b/examples/sediment_deposit/sediment_deposit.py
index 140251f5aef1a6b80d3200c04fb80ef1c11eeb96..364b3e56ea629aa09ecd56b0deb2422589d99200 100644
--- a/examples/sediment_deposit/sediment_deposit.py
+++ b/examples/sediment_deposit/sediment_deposit.py
@@ -249,7 +249,8 @@ def compute(args):
g = 1.0
Fext = SymbolicExternalForce(name='S', Fext=(0,-g*Ss),
diffusion = {S: nu_S})
- external_force = SpectralExternalForce(vorticity=vorti, dt=dt,
+ external_force = SpectralExternalForce(name='Fext',
+ 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 64238517609fb787715242becdee0a8daadaaceb..b76699dde788ead88a38d997a955466648d25576 100644
--- a/hysop/backend/device/opencl/operator/external_force.py
+++ b/hysop/backend/device/opencl/operator/external_force.py
@@ -1,19 +1,28 @@
import primefac, functools
+import sympy as sm
+
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.external_force import SpectralExternalForceOperatorBase
-from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
+from hysop.tools.types import first_not_None, to_tuple, check_instance
from hysop.core.graph.graph import op_apply
from hysop.core.memory.memory_request import MemoryRequest
+from hysop.fields.continuous_field import Field, ScalarField
+from hysop.parameters.tensor_parameter import TensorParameter
+from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.backend.device.opencl.opencl_fft import OpenClFFT
+from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
from hysop.backend.device.codegen.base.variables import dtype_to_ctype
from hysop.symbolic import local_indices_symbols
from hysop.symbolic.relational import Assignment
from hysop.symbolic.complex import ComplexMul
from hysop.symbolic.misc import Select
+from hysop.symbolic.field import AppliedSymbolicField, SymbolicField, curl, laplacian
+from hysop.symbolic.parameter import SymbolicScalarParameter
+from hysop.symbolic.spectral import AppliedSpectralTransform
+from hysop.operator.base.external_force import ExternalForce, SpectralExternalForceOperatorBase
class OpenClSpectralExternalForce(SpectralExternalForceOperatorBase, OpenClSymbolic):
"""
@@ -21,14 +30,290 @@ class OpenClSpectralExternalForce(SpectralExternalForceOperatorBase, OpenClSymbo
"""
@debug
- def __init__(self, **kwds):
- super(OpenClSpectralExternalForce, self).__init__(**kwds)
+ def __init__(self, Fext, **kwds):
+ check_instance(Fext, SymbolicExternalForce)
+ super(OpenClSpectralExternalForce, self).__init__(Fext=Fext, **kwds)
- @debug
- def setup(self, work):
- super(OpenClSpectralExternalForce, self).setup(work)
- @op_apply
- def apply(self, **kwds):
- super(OpenClSpectralExternalForce, self).apply(**kwds)
+
+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)
+ 2) Fext = (Rs*S+C)*e_y
+ Fext = (0, -Rs*S.s()+C.s())
+ """
+ Fext = to_tuple(Fext)
+ dim = len(Fext)
+
+ Fext = list(Fext)
+ for i,e in enumerate(Fext):
+ if isinstance(e, type(None)):
+ Fext[i] = 0 # curl(const) = 0
+ msg='Expression "{}" contains a SymbolicField, did you forget to apply it ?'
+ msg=msg.format(e)
+ 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(), k.short_description()
+ assert isinstance(v, ScalarParameter)
+ self._diffusion = diffusion
+
+ def input_fields(self):
+ return set(map(lambda f: f.field, self._extract_objects(AppliedSymbolicField)))
+ def output_fields(self):
+ return set(self._diffusion.keys())
+ def input_params(self):
+ p0 = set(map(lambda p: p.parameter, self._extract_objects(SymbolicScalarParameter)))
+ p1 = set(self._diffusion.values())
+ return p0.union(p1)
+ def output_params(self):
+ return set()
+
+ def initialize(self, op):
+ tg = op.new_transform_group()
+ fft_fields = tuple(self.input_fields())
+ forward_transforms = {}
+ backward_transforms = {}
+ for Si in fft_fields:
+ Fi = tg.require_forward_transform(Si)
+ if (Si in self.diffusion):
+ Bi = tg.require_backward_transform(Si)
+ else:
+ Bi = None
+ forward_transforms[Si] = Fi
+ backward_transforms[Si] = Bi
+
+ force_backward_transforms = {}
+ for Fi in op.force.fields:
+ force_backward_transforms[Fi] = tg.require_backward_transform(Fi,
+ custom_input_buffer='B0')
+
+ frame = None
+ fft_expressions = ()
+ for e in self.Fext:
+ if isinstance(e, sm.Basic):
+ efields = tuple(e.atoms(AppliedSymbolicField))
+ for sf in efields:
+ field = sf.field
+ assert field in forward_transforms, field.name
+ if field in self._diffusion:
+ assert field in backward_transforms, field.name
+ replace = {sf:forward_transforms[sf.field].s for sf in efields}
+ frame = replace.values()[0].frame
+ e = e.xreplace(replace)
+ fft_expressions += (e,)
+
+ if (frame is None):
+ msg='Could not extract frame from expressions.'
+ raise RuntimeError(frame)
+ fft_expressions = to_tuple(curl(fft_expressions, frame))
+
+ self.tg = tg
+ self.forward_transforms = forward_transforms
+ self.backward_transforms = backward_transforms
+ self.force_backward_transforms = force_backward_transforms
+ self.fft_expressions = fft_expressions
+ self.compute_required_kernels(op)
+
+ def compute_required_kernels(self, op):
+ dts = op.dt.s
+ forces = op.force.s()
+ diffusion_kernels = {}
+ for (f, nu) in self.diffusion.iteritems():
+ nus = nu.s
+ kname = 'filter_diffusion_{}d_{}'.format(f.dim, f.name)
+ Ft = self.forward_transforms[f]
+ Fs = Ft.output_symbolic_array('{}_hat'.format(f.name))
+ E = 0
+ Wn = self.tg.push_expressions(laplacian(Ft.s, Ft.s.frame))
+ for Wi in Wn:
+ Wi = self.tg._indexed_wave_numbers[Wi]
+ E += Wi
+ expr = Assignment(Fs, Fs / (1 - nus*dts*E))
+ op.require_symbolic_kernel(kname, expr)
+ diffusion_kernels[f] = kname
+
+ force_kernels = ()
+ vorticity_kernels = ()
+ names = ('dWx', 'dWy', 'dWz')
+ assert len(op.vorticity.fields)==len(op.force.fields)==len(self.fft_expressions)<=len(names)
+ for (Fi,Wi,name,e) in zip(
+ op.force.fields,
+ op.vorticity.fields,
+ names, self.fft_expressions):
+ if (e==0):
+ force_kernels += (None,)
+ vorticity_kernels += (None,)
+ continue
+
+ Fis = Fi.s()
+ Fi_hat = self.force_backward_transforms[Fi]
+ Fi_buf = Fi_hat.input_symbolic_array('{}_hat'.format(name))
+ Wn = self.tg.push_expressions(Assignment(Fi_hat, e))
+
+ msg='Could not extract transforms.'
+ try:
+ transforms = e.atoms(AppliedSpectralTransform)
+ except AttributeError:
+ raise RuntimeError(msg)
+ assert len(transforms)>=1, msg
+
+ fft_buffers = {Ft: Ft.output_symbolic_array('{}_hat'.format(Ft.field.name))}
+ wavenumbers = {Wi: self.tg._indexed_wave_numbers[Wi] for Wi in Wn}
+
+ replace = {}
+ replace.update(fft_buffers)
+ replace.update(wavenumbers)
+ expr = e.xreplace(replace)
+ expr = Assignment(Fi_buf, expr)
+
+ kname = 'compute_{}'.format(name)
+ op.require_symbolic_kernel(kname, expr)
+ force_kernels += (kname,)
+
+ Wis = Wi.s()
+ expr = Assignment(Wis, Wis + dts*Fis)
+ kname = 'update_{}'.format(Wi.var_name)
+ op.require_symbolic_kernel(kname, expr)
+ vorticity_kernels += (kname,)
+
+ assert len(diffusion_kernels) == len(self.diffusion)
+ assert len(force_kernels) == op.vorticity.nb_components == len(vorticity_kernels)
+ self.diffusion_kernel_names = diffusion_kernels
+ self.force_kernel_names = force_kernels
+ self.vorticity_kernel_names = vorticity_kernels
+
+ def discretize(self, op):
+ pass
+
+ def get_mem_requests(self, op):
+ requests = {}
+ for (Ft, Bt) in zip(self.forward_transforms.values(),
+ self.backward_transforms.values()):
+ if (Bt is not None):
+ assert (Ft.backend is Bt.backend)
+ assert (Ft.output_dtype == Bt.input_dtype), (Ft.output_dtype, Bt.input_dtype)
+ assert (Ft.output_shape == Bt.input_shape), (Ft.output_shape, Bt.input_shape)
+ shape = Ft.output_shape
+ dtype = Ft.output_dtype
+ request = MemoryRequest(backend=self.tg.backend, dtype=dtype,
+ shape=shape, nb_components=1,
+ alignment=op.min_fft_alignment)
+ name = '{}_hat'.format(Ft.field.name)
+ requests[name] = request
+ return requests
+
+ def pre_setup(self, op, work):
+ for (Ft, Bt) in zip(self.forward_transforms.values(),
+ self.backward_transforms.values()):
+ dtmp, = work.get_buffer(op, '{}_hat'.format(Ft.field.name))
+ Ft.configure_output_buffer(dtmp)
+ if (Bt is not None):
+ Bt.configure_input_buffer(dtmp)
+
+ def post_setup(self, op, work):
+ diffusion_kernels = {}
+ force_kernels = {}
+ vorticity_kernels = {}
+ ghost_exchangers = {}
+
+ queue = self.tg.backend.cl_env.default_queue
+ def build_launcher(knl, update_params):
+ def kernel_launcher(knl=knl, update_params=update_params, queue=queue):
+ kwds = update_params()
+ return knl(queue=queue, **kwds)
+ return kernel_launcher
+
+ for (field, kname) in self.diffusion_kernel_names.iteritems():
+ dfield = op.get_discrete_field(field)
+ knl, update_params = op.symbolic_kernels[kname]
+ diffusion_kernels[field] = build_launcher(knl, update_params)
+ ghost_exchangers[field] = dfield.build_ghost_exchanger(queue=queue)
+
+ for i, (kname0, kname1) in enumerate(zip(
+ self.force_kernel_names, self.vorticity_kernel_names)):
+ if (kname0 is None):
+ assert (kname1 is None)
+ continue
+ Wi = op.vorticity[i]
+ Fi = op.force[i]
+ dWi = op.dW[i]
+
+ knl, update_params = op.symbolic_kernels[kname0]
+ force_kernels[(Fi,Wi)] = build_launcher(knl, update_params)
+
+ knl, update_params = op.symbolic_kernels[kname1]
+ vorticity_kernels[(Fi,Wi)] = build_launcher(knl, update_params)
+
+ ghost_exchangers[Wi] = dWi.build_ghost_exchanger(queue=queue)
+
+ def compute_statistics():
+ pass
+ self.compute_statistics = compute_statistics
+
+ def apply(self, op, **kwds):
+ for (field, Ft) in self.forward_transforms.iteritems():
+ evt = Ft()
+ if (field in self.backward_transforms):
+ evt = self.diffusion_kernels[field]()
+ evt = self.backward_transforms[field]()
+ evt = self.ghost_exchangers[field]()
+
+ for (Fi,Wi) in self.force_kernels.keys():
+ evt = self.force_kernel[Wi]()
+ evt = self.force_backward_transforms[Fi]()
+ if op.compute_statistics:
+ evt = self.compute_statistics()
+ evt = self.vorticity_kernels[Wi]()
+ evt = self.ghost_exchangers[Wi]()
+
+ def _extract_objects(self, obj_type):
+ objs = set()
+ for e in self.Fext:
+ 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
+
diff --git a/hysop/core/graph/computational_operator.py b/hysop/core/graph/computational_operator.py
index 5c3e192bfb7612fc6d6e8284b194e88b346ac01a..3adf180d71a44e4969a73475e17408f1e8b78d06 100644
--- a/hysop/core/graph/computational_operator.py
+++ b/hysop/core/graph/computational_operator.py
@@ -515,7 +515,10 @@ class ComputationalGraphOperator(ComputationalGraphNode):
Automatically honour temporary field memory requests.
"""
+ self.allocate_tmp_fields(work)
super(ComputationalGraphOperator, self).setup(work)
+
+ def allocate_tmp_fields(self, work):
for dfield in self.discrete_fields:
if dfield.is_tmp:
req_id = 'tmp_{}_{}'.format(dfield.name, dfield.tag)
diff --git a/hysop/operator/base/external_force.py b/hysop/operator/base/external_force.py
index 0fd15eaee75424528f6ce35d7bb06fa7b86b5115..f33c0d4b2b1868911cfaf09cd33e3a4c274c115a 100644
--- a/hysop/operator/base/external_force.py
+++ b/hysop/operator/base/external_force.py
@@ -1,21 +1,79 @@
-
-
-from abc import abstractmethod
+from abc import ABCMeta, abstractmethod
from hysop.constants import SpectralTransformAction
from hysop.tools.numpywrappers import npw
from hysop.tools.types import check_instance, first_not_None, to_tuple
from hysop.tools.decorators import debug
from hysop.tools.numerics import float_to_complex_dtype
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import Field, ScalarField
from hysop.operator.base.spectral_operator import SpectralOperatorBase
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.fields.continuous_field import Field
from hysop.symbolic.relational import Assignment
from hysop.core.memory.memory_request import MemoryRequest
-from hysop.operator.external_force import ExternalForce
from hysop.parameters.tensor_parameter import TensorParameter
from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.tools.interface import NamedObjectI
+
+class ExternalForce(NamedObjectI):
+ """Interface to implement a custom external force."""
+ __metaclass__ = ABCMeta
+
+ def __init__(self, name, dim, Fext, **kwds):
+ super(ExternalForce, self).__init__(name=name, **kwds)
+
+ check_instance(dim, int)
+ self._dim = dim
+ self._Fext = Fext
+
+ for f in self.input_fields():
+ msg='Expected field dimension to be {} but got {} from field {}.'
+ msg=msg.format(dim, f.dim, f.short_description())
+ assert f.dim == dim, msg
+
+
+ @property
+ def Fext(self):
+ return self._Fext
+ @property
+ def dim(self):
+ return self._dim
+ @property
+ def diffusion(self):
+ return self._diffusion
+
+ @abstractmethod
+ def input_fields(self):
+ pass
+ @abstractmethod
+ def output_fields(self):
+ pass
+ @abstractmethod
+ def input_params(self):
+ pass
+ @abstractmethod
+ def output_params(self):
+ pass
+
+ @abstractmethod
+ def initialize(self, op):
+ pass
+ @abstractmethod
+ def discretize(self, op):
+ pass
+ @abstractmethod
+ def get_mem_requests(self, op):
+ pass
+ @abstractmethod
+ def pre_setup(self, op, work):
+ pass
+ @abstractmethod
+ def post_setup(self, op, work):
+ pass
+ @abstractmethod
+ def apply(self, op, **kwds):
+ pass
+
class SpectralExternalForceOperatorBase(SpectralOperatorBase):
"""
@@ -34,11 +92,11 @@ class SpectralExternalForceOperatorBase(SpectralOperatorBase):
vorticity: 2 1 | 3 3
What is computed:
- tmp = curl(Fext) by using a spectral backend
- Fmin = min(tmp)
- Fmax = max(tmp)
+ force = curl(Fext) by using a spectral backend
+ Fmin = min(force)
+ Fmax = max(force)
Finf = max(abs(Fmin), abs(Fmax))
- W += dt*tmp
+ W += dt*force
where Fext is computed from user given ExternalForce.
@@ -115,14 +173,20 @@ class SpectralExternalForceOperatorBase(SpectralOperatorBase):
assert Fmin.shape==pshape, msg.format(Finf.shape, 'Finf')
Finf = Finf.view((0,)) if (dim==2) else Finf
+ compute_statistics = (Fmin is not None)
+ compute_statistics |= (Fmax is not None)
+ compute_statistics |= (Finf is not None)
+
# input and output fields
input_fields = Fext.input_fields()
+ check_instance(input_fields, set, values=ScalarField)
input_fields.add(vorticity)
output_fields = Fext.output_fields()
+ check_instance(output_fields, set, values=ScalarField)
output_fields.add(vorticity)
- input_params = Fext.input_params()
+ input_params = Fext.input_params()
input_params.add(dt)
output_params = Fext.output_params()
@@ -133,6 +197,12 @@ class SpectralExternalForceOperatorBase(SpectralOperatorBase):
output_fields = {f: self.get_topo_descriptor(variables, f) for f in output_fields}
input_params = {p.name: p for p in input_params}
output_params = {p.name: p for p in output_params}
+
+ # TODO share tmp buffers for the whole tensor
+ force = vorticity.tmp_like(name='Fext', ghosts=0)
+ for (Fi, Wi) in zip(force.fields, vorticity.fields):
+ input_fields[Fi] = self.get_topo_descriptor(variables, Wi)
+ output_fields[Fi] = self.get_topo_descriptor(variables, Wi)
super(SpectralExternalForceOperatorBase, self).__init__(
input_fields=input_fields, output_fields=output_fields,
@@ -140,18 +210,47 @@ class SpectralExternalForceOperatorBase(SpectralOperatorBase):
**kwds)
self.vorticity = vorticity
- self.Fext = Fext
+ self.Fext = Fext
+ self.force = force
+ self.dt = dt
+
self.dim = dim
self.w_components = w_components
self.f_components = f_components
+
self.Fmin = Fmin
self.Fmax = Fmax
self.Finf = Finf
+ self.compute_statistics = compute_statistics
+
+###################
+# from now on, we delegate everything to the ExternalForce implementation
+###################
+ def initialize(self, **kwds):
+ if self.initialized:
+ return
+ self.Fext.initialize(self)
+ return super(SpectralExternalForceOperatorBase, self).initialize(**kwds)
@debug
def discretize(self):
if self.discretized:
return
super(SpectralExternalForceOperatorBase, self).discretize()
- self.dW = self.get_input_discrete_field(self.vorticity)
+ self.dW = self.get_input_discrete_field(self.vorticity)
+ self.Fext.discretize(self)
+
+ def get_work_properties(self):
+ requests = super(SpectralExternalForceOperatorBase, self).get_work_properties()
+ for (name, request) in self.Fext.get_mem_requests(self).iteritems():
+ requests.push_mem_request(name, request)
+ return requests
+
+ def setup(self, work):
+ self.Fext.pre_setup(self, work)
+ super(SpectralExternalForceOperatorBase, self).setup(work)
+ self.Fext.post_setup(self, work)
+
+ def apply(self, **kwds):
+ self.Fext.apply(self, **kwds)
diff --git a/hysop/operator/base/spectral_operator.py b/hysop/operator/base/spectral_operator.py
index 4d7a1a666b68e51d5608c902c6a353d2036adffd..ba5a2782877075878ea2d0c5e38adfdebf3c7a5b 100644
--- a/hysop/operator/base/spectral_operator.py
+++ b/hysop/operator/base/spectral_operator.py
@@ -220,6 +220,7 @@ class SpectralOperatorBase(object):
def setup(self, work):
+ self.allocate_tmp_fields(work)
for tg in self.transform_groups.values():
tg.setup(work=work)
super(SpectralOperatorBase, self).setup(work=work)
@@ -1125,12 +1126,6 @@ class PlannedSpectralTransform(object):
self._backend = dfield.backend
- # bind field buffer to input or output
- if is_forward:
- self.configure_input_buffer(dfield.sbuffer[dfield.compute_slices])
- else:
- self.configure_output_buffer(dfield.sbuffer[dfield.compute_slices])
-
if self.DEBUG:
def axis_format(info):
prefix='\n'+' '*4
@@ -1526,6 +1521,13 @@ class PlannedSpectralTransform(object):
TMP, = work.get_buffer(op, TMP_tag, handle=True)
except ValueError:
TMP = None
+
+ # bind field buffer to input or output
+ dfield = self.dfield
+ if is_forward:
+ self.configure_input_buffer(dfield.sbuffer[dfield.compute_slices])
+ else:
+ self.configure_output_buffer(dfield.sbuffer[dfield.compute_slices])
# bind group buffer to input or output if required.
custom_input_buffer = self._custom_input_buffer
diff --git a/hysop/operator/external_force.py b/hysop/operator/external_force.py
index bb77ad21c5bfa09be1ba8c136a8ae342c36d7c51..0709bbfbe17bb743a2cce46336f34a92358116b4 100644
--- a/hysop/operator/external_force.py
+++ b/hysop/operator/external_force.py
@@ -1,10 +1,6 @@
-import sympy as sm
-from abc import ABCMeta, abstractmethod
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.fields.continuous_field import Field, ScalarField
from hysop.parameters.tensor_parameter import TensorParameter
from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.operator.min_max import MinMaxFieldStatisticsBase
@@ -12,131 +8,8 @@ from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeF
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
-
- def __init__(self, name, dim, Fext, **kwds):
- super(ExternalForce, self).__init__(name=name, **kwds)
-
- check_instance(dim, int)
- self._dim = dim
- self._Fext = Fext
-
- for f in self.input_fields():
- msg='Expected field dimension to be {} but got {} from field {}.'
- msg=msg.format(dim, f.dim, f.short_description())
- assert f.dim == dim, msg
-
-
- @property
- def Fext(self):
- return self._Fext
- @property
- def dim(self):
- return self._dim
- @property
- def diffusion(self):
- return self._diffusion
-
- @abstractmethod
- def input_fields(self):
- pass
- @abstractmethod
- def output_fields(self):
- pass
- @abstractmethod
- def input_params(self):
- pass
- @abstractmethod
- def output_params(self):
- pass
-
-
-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)
- 2) Fext = (Rs*S+C)*e_y
- Fext = (0, -Rs*S.s()+C.s())
- """
- Fext = to_tuple(Fext)
- dim = len(Fext)
-
- Fext = list(Fext)
- for i,e in enumerate(Fext):
- if isinstance(e, (type(None), int, long, float)):
- Fext[i] = 0 # curl(const) = 0
- msg='Expression "{}" contains a SymbolicField, did you forget to apply it ?'
- msg=msg.format(e)
- 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(), k.short_description()
- assert isinstance(v, ScalarParameter)
- self._diffusion = diffusion
-
- def input_fields(self):
- return set(map(lambda f: f.field, self._extract_objects(AppliedSymbolicField)))
- def output_fields(self):
- return set(self._diffusion.keys())
- def input_params(self):
- p0 = set(map(lambda p: p.parameter, self._extract_objects(SymbolicScalarParameter)))
- p1 = set(self._diffusion.values())
- return p0.union(p1)
- def output_params(self):
- return set()
-
- def _extract_objects(self, obj_type):
- objs = set()
- for e in self.Fext:
- 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):
"""
@@ -217,6 +90,7 @@ class SpectralExternalForce(ComputationalGraphNodeFrontend):
Fext has 3 components
Expected parameters are TensorParameter of shape (3,)
"""
+ from hysop.operator.base.external_force import ExternalForce
base_kwds = first_not_None(base_kwds, {})
check_instance(vorticity, Field)
check_instance(Fext, ExternalForce)
diff --git a/hysop/operators.py b/hysop/operators.py
index ba66000fe6664c7af26e9a90903dd32859a1db8a..a68cb3fe02ae2a862f0e03abd348e451f6b5eac4 100644
--- a/hysop/operators.py
+++ b/hysop/operators.py
@@ -32,7 +32,8 @@ from hysop.operator.min_max import MinMaxFieldStatistics,
from hysop.operator.gradient import Gradient, MinMaxGradientStatistics
from hysop.operator.curl import Curl, SpectralCurl
-from hysop.operator.external_force import SymbolicExternalForce, SpectralExternalForce
+from hysop.operator.external_force import SpectralExternalForce
+from hysop.backend.device.opencl.operator.external_force import SymbolicExternalForce
from hysop.numerics.splitting.strang import StrangSplitting
from hysop.operator.directional.symbolic_dir import DirectionalSymbolic