From 2b0f162d1390adbec9013a3f3fe0aeec2b929b5a Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@univ-pau.fr>
Date: Wed, 27 Jun 2018 15:10:07 +0200
Subject: [PATCH] Add python backend for enstrophy operator
---
.../device/opencl/operator/enstrophy.py | 15 +++-----
.../backend/host/python/operator/enstrophy.py | 37 ++++++++++++++++++
hysop/operator/base/enstrophy.py | 31 +++++++--------
hysop/operator/enstrophy.py | 38 +++++++++----------
4 files changed, 78 insertions(+), 43 deletions(-)
create mode 100644 hysop/backend/host/python/operator/enstrophy.py
diff --git a/hysop/backend/device/opencl/operator/enstrophy.py b/hysop/backend/device/opencl/operator/enstrophy.py
index 0a8f8ec3c..766bdb3c5 100644
--- a/hysop/backend/device/opencl/operator/enstrophy.py
+++ b/hysop/backend/device/opencl/operator/enstrophy.py
@@ -1,4 +1,3 @@
-
from hysop.deps import sm
from hysop.constants import DirectionLabels
from hysop.backend.device.opencl.opencl_array_backend import OpenClArrayBackend
@@ -23,13 +22,13 @@ class OpenClEnstrophy(EnstrophyBase, OpenClSymbolic):
rho = self.rho
rhos = (rho.s if (rho is not None) else sm.Integer(1))
-
+
Ws = self.vorticity.s()
WdotW = self.WdotW.s()
expr = Assignment(WdotW, rhos*npw.dot(Ws, Ws))
self.require_symbolic_kernel('WdotW', expr)
-
+
@debug
def get_field_requirements(self):
# force 0 ghosts for the reduction (pyopencl reduction kernel)
@@ -38,21 +37,19 @@ class OpenClEnstrophy(EnstrophyBase, OpenClSymbolic):
if field is self.WdotW:
req.max_ghosts = (0,)*self.WdotW.dim
return requirements
-
+
@debug
def discretize(self):
super(OpenClEnstrophy, self).discretize()
assert (self.dWdotW.ghosts == 0).all()
- self.coeff = npw.prod(self.dWdotW.space_step)
- self.coeff /= (self.rho_0 * npw.prod(self.dWdotW.domain.length))
-
+
@debug
def setup(self, work):
super(OpenClEnstrophy, self).setup(work)
(self.WdotW_kernel, self.WdotW_update_parameters) = self.symbolic_kernels['WdotW']
- self.sum_kernel = self.dWdotW.backend.nansum(a=self.dWdotW[0],
+ self.sum_kernel = self.dWdotW.backend.nansum(a=self.dWdotW[0],
build_kernel_launcher=True, async=True)
-
+
@op_apply
def apply(self, **kwds):
queue = self.cl_env.default_queue
diff --git a/hysop/backend/host/python/operator/enstrophy.py b/hysop/backend/host/python/operator/enstrophy.py
new file mode 100644
index 000000000..3ec1080dd
--- /dev/null
+++ b/hysop/backend/host/python/operator/enstrophy.py
@@ -0,0 +1,37 @@
+"""
+@file enstrophy.py
+Enstrophy solver python backend.
+"""
+from hysop.tools.decorators import debug
+from hysop.core.graph.graph import op_apply
+from hysop.operator.base.enstrophy import EnstrophyBase
+from hysop.backend.host.host_operator import HostOperator
+from hysop.tools.numpywrappers import npw
+
+
+class PythonEnstrophy(EnstrophyBase, HostOperator):
+ """Compute enstrophy of the given vorticity field."""
+ @debug
+ def __init__(self, **kwds):
+ super(PythonEnstrophy, self).__init__(**kwds)
+ assert self.mpi_params.size == 1
+
+ @op_apply
+ def apply(self, **kwds):
+ print self.dWin.data[0].shape, self.dWdotW.data[0].shape
+ print len(self.dWin.data), len(self.dWdotW.data)
+
+ # Compute W dot W
+ wd = self.dWin
+ wdotw = self.dWdotW[0]
+ nbc = wd.nb_components
+ ic = wd.compute_slices
+ wdotw[...] = 0.
+ for i in xrange(nbc):
+ wdotw[...] += (wd[i][ic] ** 2)
+
+ # Compute enstrophy
+ local_enstrophy = self.coeff * npw.real_sum(wdotw)
+ # TODO reduce over all mpi process
+
+ self.enstrophy.value = local_enstrophy
diff --git a/hysop/operator/base/enstrophy.py b/hysop/operator/base/enstrophy.py
index ea1e81478..6f46e2750 100644
--- a/hysop/operator/base/enstrophy.py
+++ b/hysop/operator/base/enstrophy.py
@@ -1,10 +1,9 @@
-
-
from abc import ABCMeta
from hysop.tools.types import check_instance, to_tuple, first_not_None
from hysop.tools.decorators import debug
from hysop.fields.continuous_field import Field
+from hysop.tools.numpywrappers import npw
from hysop.core.memory.memory_request import MemoryRequest
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
@@ -22,13 +21,13 @@ class EnstrophyBase(object):
variables, name=None, pretty_name=None, **kwds):
"""
Initialize a enstrophy operator operating on CartesianTopologyDescriptors.
-
+
vorticity: Field
Input continuous vorticity field.
enstrophy: ScalarParameter
Enstrophy scalar output parameter.
rho: Field, optional
- Input continuous density field, if not given,
+ Input continuous density field, if not given,
defaults to 1.0 on the whole domain.
rho_0: float, optional
Reference density, defaults to 1.0.
@@ -36,40 +35,40 @@ class EnstrophyBase(object):
Output continuous field, will contain W.W (term before integration).
If WdotW is given, WdotW will contain W.W else this will be
a temporary field only usable during this operator's apply method.
- Should have nb_components=1 and the same domain and discretization
+ Should have nb_components=1 and the same domain and discretization
as the vorticity.
variables: dict
dictionary of fields as keys and topologies as values.
kwds:
- Keywords arguments that will be passed towards implementation
+ Keywords arguments that will be passed towards implementation
poisson operator __init__.
"""
-
+
check_instance(vorticity, Field)
check_instance(enstrophy, ScalarParameter)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(WdotW, Field, allow_none=True)
check_instance(rho, Field, allow_none=True)
check_instance(rho_0, float)
-
+
W = vorticity.pretty_name.decode('utf-8')
wdotw = u'{}\u22c5{}'.format(W,W)
if (WdotW is None):
assert (vorticity in variables), variables
WdotW = vorticity.tmp_like(name='WdotW', pretty_name=wdotw, nb_components=1)
variables[WdotW] = variables[vorticity]
-
+
input_fields = { vorticity: variables[vorticity] }
output_fields = { WdotW: variables[WdotW] }
output_params = { enstrophy.name: enstrophy }
if (rho is not None):
input_fields[rho] = variables[rho]
-
+
default_name = 'enstrophy'
default_pname = u'\u222b{}'.format(wdotw).encode('utf-8')
-
+
pretty_name = first_not_None(pretty_name, name, default_pname)
name = first_not_None(name, default_name)
@@ -78,11 +77,11 @@ class EnstrophyBase(object):
self.rho_0 = rho_0
self.WdotW = WdotW
self.enstrophy = enstrophy
-
- super(EnstrophyBase, self).__init__(input_fields=input_fields,
- output_fields=output_fields, output_params=output_params,
+
+ super(EnstrophyBase, self).__init__(input_fields=input_fields,
+ output_fields=output_fields, output_params=output_params,
name=name, pretty_name=pretty_name, **kwds)
-
+
@debug
def discretize(self):
if self.discretized:
@@ -91,3 +90,5 @@ class EnstrophyBase(object):
self.dWin = self.input_discrete_fields[self.vorticity]
self.dWdotW = self.output_discrete_fields[self.WdotW]
+ self.coeff = npw.prod(self.dWdotW.space_step)
+ self.coeff /= (self.rho_0 * npw.prod(self.dWdotW.domain.length))
diff --git a/hysop/operator/enstrophy.py b/hysop/operator/enstrophy.py
index a8637dcd4..09b2a3276 100644
--- a/hysop/operator/enstrophy.py
+++ b/hysop/operator/enstrophy.py
@@ -1,5 +1,3 @@
-
-
"""
@file enstrophy.py
Enstrophy solver frontend.
@@ -17,37 +15,40 @@ from hysop.parameters.scalar_parameter import ScalarParameter
class Enstrophy(ComputationalGraphNodeFrontend):
"""
Interface computing enstrophy.
- Available implementations are:
- *OPENCL (gpu based implementation)
+ Available implementations are:
+ *OPENCL (gpu based implementation) (default)
+ *PYTHON
"""
-
+
@classmethod
def implementations(cls):
from hysop.backend.device.opencl.operator.enstrophy import OpenClEnstrophy
+ from hysop.backend.host.python.operator.enstrophy import PythonEnstrophy
implementations = {
- Implementation.OPENCL: OpenClEnstrophy
+ Implementation.OPENCL: OpenClEnstrophy,
+ Implementation.PYTHON: PythonEnstrophy
}
return implementations
-
+
@classmethod
def default_implementation(cls):
return Implementation.OPENCL
-
+
@debug
- def __init__(self, vorticity, enstrophy, variables,
- rho=None, rho_0=1.0, WdotW=None,
+ def __init__(self, vorticity, enstrophy, variables,
+ rho=None, rho_0=1.0, WdotW=None,
implementation=None, base_kwds=None, **kwds):
"""
Initialize a Enstrophy operator frontend.
- Enstrophy is the scaled volume average of rho*(W.W) on the domain where . represents
+ Enstrophy is the scaled volume average of rho*(W.W) on the domain where . represents
the vector dot product).
in: W (vorticity field)
rho (density field, optional, defaults to 1.0 everywhere)
out: E = 1.0/(2*V*rho_0) * integral(rho*(W.W)) => enstrophy (scalar parameter)
where V is the domain volume, rho_0 the reference density.
-
+
Parameters
----------
vorticity: Field
@@ -55,7 +56,7 @@ class Enstrophy(ComputationalGraphNodeFrontend):
enstrophy: ScalarParameter
Enstrophy scalar output parameter.
rho: Field, optional
- Input continuous density field, if not given,
+ Input continuous density field, if not given,
defaults to 1.0 on the whole domain.
rho_0: float, optional
Reference density, defaults to 1.0.
@@ -73,12 +74,12 @@ class Enstrophy(ComputationalGraphNodeFrontend):
Base class keywords arguments.
If None, an empty dict will be passed.
kwds:
- Extra keywords arguments that will be passed towards implementation
+ Extra keywords arguments that will be passed towards implementation
enstrophy operator __init__.
Notes
-----
- An Enstrophy operator implementation should at least support
+ An Enstrophy operator implementation should at least support
the hysop.operator.base.enstrophy.EnstrophyBase interface.
"""
base_kwds = base_kwds or dict()
@@ -90,9 +91,8 @@ class Enstrophy(ComputationalGraphNodeFrontend):
check_instance(WdotW, Field, allow_none=True)
check_instance(rho, Field, allow_none=True)
check_instance(rho_0, float)
-
- super(Enstrophy, self).__init__(vorticity=vorticity, rho=rho,
+
+ super(Enstrophy, self).__init__(vorticity=vorticity, rho=rho,
enstrophy=enstrophy, WdotW=WdotW, rho_0=rho_0,
- variables=variables, base_kwds=base_kwds,
+ variables=variables, base_kwds=base_kwds,
implementation=implementation, **kwds)
-
--
GitLab