From 68a824db7c426d716f0dfa38de4ff523b0256178 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 31 Mar 2019 12:25:16 +0200
Subject: [PATCH] fixed analytic, fixed levelset example
---
examples/scalar_advection/levelset.py | 89 ++++++++++++++-----
.../device/opencl/operator/analytic.py | 2 +-
.../device/opencl/operator/integrate.py | 3 +-
.../backend/host/python/operator/analytic.py | 16 ++--
hysop/operator/base/enstrophy.py | 1 +
hysop/operator/base/integrate.py | 3 +-
hysop/operator/integrate.py | 5 +-
7 files changed, 87 insertions(+), 32 deletions(-)
diff --git a/examples/scalar_advection/levelset.py b/examples/scalar_advection/levelset.py
index 9246d2060..21f3ebec7 100644
--- a/examples/scalar_advection/levelset.py
+++ b/examples/scalar_advection/levelset.py
@@ -3,11 +3,11 @@ import sympy as sm
pi = np.pi
def compute(args):
- from hysop import Field, Box, Simulation, Problem, \
+ from hysop import Field, Box, Simulation, Problem, IOParams, \
ScalarParameter, MPIParams, CartesianDiscretization, CartesianTopology
from hysop.constants import Implementation, Backend
from hysop.operators import DirectionalAdvection, StrangSplitting, Integrate, \
- AnalyticField, Advection
+ AnalyticField, Advection, HDF_Writer
from hysop.methods import Remesh, TimeIntegrator, ComputeGranularity, \
Interpolation, StrangOrder
import hysop.numerics.odesolvers.runge_kutta as rk
@@ -51,8 +51,9 @@ def compute(args):
device_id=args.cl_device_id)
extra_op_kwds['cl_env'] = cl_env
method[OpenClKernelConfig] = args.opencl_kernel_config
+ backend = Backend.OPENCL
elif (impl in (Implementation.PYTHON, Implementation.FORTRAN)):
- pass
+ backend = Backend.HOST
else:
msg='Unknown implementation \'{}\'.'.format(impl)
raise ValueError(msg)
@@ -63,6 +64,11 @@ def compute(args):
scalar = Field(domain=box, name='S', nb_components=1, dtype=args.dtype)
vol = ScalarParameter('volume', dtype=args.dtype)
+ if (snpts == npts):
+ mask = velo[1]
+ else:
+ mask = scalar.field_like('mask')
+
# Setup operator method dictionnary
method.update({
TimeIntegrator: args.time_integrator,
@@ -90,9 +96,15 @@ def compute(args):
data[1][...] = -sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)*cos(t()*pi/3.)
(x,y,z) = coords[2]
data[2][...] = -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t()*pi/3.)
- formula = compute_velocity
+ def compute_norm(data, coords):
+ data[0][...] = np.sqrt(data[0]**2 + data[1]**2 + data[2]**2)
+ def compute_volume(data, coords, S):
+ data[0][...] = (S[0] <= 0.5)
elif (impl is Implementation.OPENCL):
+ from hysop.symbolic.relational import LogicalLE
sin, cos = sm.sin, sm.cos
+ Vs = velo.s()
+ Ss = scalar.s()
if dim == 3:
x, y, z = box.frame.coords
formula = (2.*sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)*cos(t.s*pi/3.),
@@ -102,24 +114,13 @@ def compute(args):
x, y = box.frame.coords
formula = (-sin(x*pi)**2*sin(y*pi*2)*cos(t.s*pi/3.),
sin(y*pi)**2*sin(x*pi*2)*cos(t.s*pi/3.))
+ compute_velocity = formula
+ compute_norm = (sm.sqrt(np.dot(Vs, Vs)), None, None)
+ compute_volume = LogicalLE(Ss, 0.5)
else:
msg='Unknown implementation \'{}\'.'.format(impl)
raise ValueError(msg)
-
- analytic = AnalyticField(field=velo, formula=formula,
- variables={velo: npts},
- extra_input_kwds={'t': simu.t},
- implementation=Implementation.PYTHON if (impl is Implementation.FORTRAN) else impl,
- **extra_op_kwds)
-
- volume = Integrate(name='volume',
- field=scalar,
- variables={scalar: snpts},
- parameter=vol,
- implementation=impl,
- **extra_op_kwds)
-
if (impl is Implementation.FORTRAN):
advec = Advection(name='advec',
velocity = velo,
@@ -140,12 +141,58 @@ def compute(args):
advec = StrangSplitting(splitting_dim=dim,
order=args.strang_order)
advec.push_operators(advec_dir)
+
+
+ compute_velocity = AnalyticField(name='compute_velocity', pretty_name='V(t)',
+ field=velo, formula=compute_velocity,
+ variables={velo: npts},
+ extra_input_kwds={'t': simu.t},
+ implementation=Implementation.PYTHON if (impl is Implementation.FORTRAN) else impl,
+ **extra_op_kwds)
+
+ compute_norm = AnalyticField(name='compute_norm', pretty_name='||V||',
+ field=velo, formula=compute_norm,
+ variables={velo: npts},
+ implementation=Implementation.PYTHON if (impl is Implementation.FORTRAN) else impl,
+ **extra_op_kwds)
+
+ compute_volume = AnalyticField(name='compute_volume', pretty_name='S<0.5',
+ field=mask, formula=compute_volume,
+ variables={mask: snpts, scalar: snpts},
+ extra_input_kwds={'S': scalar},
+ implementation=Implementation.PYTHON if (impl is Implementation.FORTRAN) else impl,
+ **extra_op_kwds)
+
+ # compute the volume (where S<=0.5), use Vy as tmp buffer when possible (ie. npts == snpts)
+ volume = Integrate(name='volume',
+ field=mask,
+ variables={mask: snpts},
+ parameter=vol,
+ implementation=impl,
+ allow_ghosts=True,
+ **extra_op_kwds)
+
+ # dump scalar and velocity norm
+ if (npts==snpts):
+ io_params = IOParams(filename='fields', frequency=args.dump_freq)
+ dump_fields = HDF_Writer(name='fields', io_params=io_params,
+ force_backend=backend, variables={velo[0]: snpts, scalar: snpts}, **extra_op_kwds)
+ dumpers = (dump_fields,)
+ else:
+ io_params = IOParams(filename='scalar', frequency=args.dump_freq)
+ dump_scalar = HDF_Writer(name='scalar', io_params=io_params,
+ force_backend=backend, variables={scalar: snpts}, **extra_op_kwds)
+ io_params = IOParams(filename='Vnorm', frequency=args.dump_freq)
+ dump_velocity_norm = HDF_Writer(name='Vnorm', io_params=io_params,
+ force_backend=backend, variables={velo[0]: npts}, **extra_op_kwds)
+ dumpers = (dump_scalar, dump_velocity_norm)
problem = Problem(method=method)
- problem.insert(analytic, advec, volume)
- problem.dump_inputs(fields=scalar, filename='S', frequency=args.dump_freq)
- problem.build()
+ problem.insert(compute_velocity, advec)
+ problem.insert(compute_norm, *dumpers)
+ problem.insert(compute_volume, volume)
+ problem.build(args=args)
# If a visu_rank was provided, and show_graph was set,
# display the graph on the given process rank.
diff --git a/hysop/backend/device/opencl/operator/analytic.py b/hysop/backend/device/opencl/operator/analytic.py
index 7ea4fce4e..d490f3654 100644
--- a/hysop/backend/device/opencl/operator/analytic.py
+++ b/hysop/backend/device/opencl/operator/analytic.py
@@ -34,7 +34,7 @@ class OpenClAnalyticField(OpenClCustomSymbolicOperator):
"""
formula = to_tuple(formula)
check_instance(field, Field)
- check_instance(formula, tuple, values=sm.Basic, size=field.nb_components)
+ check_instance(formula, tuple, values=(type(None),sm.Basic), size=field.nb_components)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
exprs = ()
diff --git a/hysop/backend/device/opencl/operator/integrate.py b/hysop/backend/device/opencl/operator/integrate.py
index d3feef93f..6aa31e79d 100644
--- a/hysop/backend/device/opencl/operator/integrate.py
+++ b/hysop/backend/device/opencl/operator/integrate.py
@@ -13,7 +13,8 @@ class OpenClIntegrate(IntegrateBase, OpenClOperator):
# force 0 ghosts for the reduction (pyopencl reduction kernel)
requirements = super(OpenClIntegrate, self).get_field_requirements()
topo, req = requirements.get_input_requirement(self.field)
- req.max_ghosts = (0,)*self.field.dim
+ if not self.allow_ghosts:
+ req.max_ghosts = (0,)*self.field.dim
return requirements
@debug
diff --git a/hysop/backend/host/python/operator/analytic.py b/hysop/backend/host/python/operator/analytic.py
index 4f9ffb5a2..e2969cfe0 100644
--- a/hysop/backend/host/python/operator/analytic.py
+++ b/hysop/backend/host/python/operator/analytic.py
@@ -48,15 +48,17 @@ class PythonAnalyticField(HostOperator):
assert callable(formula), type(formula)
check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(extra_input_kwds, dict, keys=str)
-
+
input_fields = {}
output_fields = { field: variables[field] }
input_params = {}
extra_kwds = {}
+ map_fields = {}
for (k,v) in extra_input_kwds.iteritems():
if isinstance(v, Field):
- input_fields[k] = v
+ input_fields[v] = variables[v]
+ map_fields[v] = k
elif isinstance(v, Parameter):
input_params[k] = v
extra_kwds[k] = v
@@ -69,6 +71,7 @@ class PythonAnalyticField(HostOperator):
self.field = field
self.formula = formula
self.extra_kwds = extra_kwds
+ self.map_fields = map_fields
@debug
def discretize(self):
@@ -77,13 +80,14 @@ class PythonAnalyticField(HostOperator):
super(PythonAnalyticField, self).discretize()
dfield = self.get_output_discrete_field(self.field)
extra_kwds = self.extra_kwds
+ map_fields = self.map_fields
assert 'data' not in extra_kwds
assert 'coords' not in extra_kwds
- extra_kwds['data'] = dfield.data
- extra_kwds['coords'] = dfield.get_attributes('mesh', 'local_mesh_coords')
- for (field, dfield) in self.input_discrete_fields:
+ extra_kwds['data'] = dfield.compute_data
+ extra_kwds['coords'] = dfield.get_attributes('compute_mesh_coords')
+ for (field, dfield) in self.input_discrete_fields.iteritems():
assert field.name not in extra_kwds, field.name
- extra_kwds[field.name] = dfield.data
+ extra_kwds[map_fields[field]] = dfield.compute_data
self.dfield = dfield
@op_apply
diff --git a/hysop/operator/base/enstrophy.py b/hysop/operator/base/enstrophy.py
index 8dd2a558f..1305461ca 100644
--- a/hysop/operator/base/enstrophy.py
+++ b/hysop/operator/base/enstrophy.py
@@ -112,3 +112,4 @@ class EnstrophyBase(object):
@classmethod
def supports_mpi(cls):
return True
+
diff --git a/hysop/operator/base/integrate.py b/hysop/operator/base/integrate.py
index d6c152323..d3b1b9396 100644
--- a/hysop/operator/base/integrate.py
+++ b/hysop/operator/base/integrate.py
@@ -17,7 +17,7 @@ class IntegrateBase(object):
__metaclass__ = ABCMeta
@debug
- def __init__(self, field, variables,
+ def __init__(self, field, variables, allow_ghosts=False,
name=None, pretty_name=None, cst=1,
parameter=None, scaling=None, **kwds):
"""
@@ -89,6 +89,7 @@ class IntegrateBase(object):
self.scaling = scaling
self.cst = cst
self.scaling_coeff = None
+ self.allow_ghosts = allow_ghosts
super(IntegrateBase, self).__init__(name=name, pretty_name=pretty_name,
input_fields=input_fields, output_params=output_params, **kwds)
diff --git a/hysop/operator/integrate.py b/hysop/operator/integrate.py
index 4c0a8e646..34de06ea2 100644
--- a/hysop/operator/integrate.py
+++ b/hysop/operator/integrate.py
@@ -36,7 +36,7 @@ class Integrate(ComputationalGraphNodeFrontend):
@debug
def __init__(self, field, variables,
- parameter=None, scaling=None,
+ parameter=None, scaling=None,
implementation=None, base_kwds=None, **kwds):
"""
Initialize a Integrate operator frontend.
@@ -97,4 +97,5 @@ class Integrate(ComputationalGraphNodeFrontend):
raise RuntimeError(msg)
super(Integrate, self).__init__(field=field, variables=variables,
- parameter=parameter, scaling=scaling, base_kwds=base_kwds, **kwds)
+ parameter=parameter, scaling=scaling,
+ base_kwds=base_kwds, **kwds)
--
GitLab