From a2693b7ebad3c6ed1858c8814cb4cdaa616cfae4 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 3 Feb 2019 19:36:59 +0100
Subject: [PATCH] working python example
---
examples/multiresolution/scalar_advection.py | 37 +++--
.../host/python/operator/spatial_filtering.py | 27 ++++
hysop/core/graph/computational_node.py | 2 +-
hysop/core/graph/computational_operator.py | 19 ++-
hysop/core/graph/graph_builder.py | 32 +++--
hysop/numerics/remesh/kernel_generator.py | 24 +++-
hysop/numerics/remesh/remesh.py | 2 +-
hysop/operator/base/spatial_filtering.py | 129 +++++++++++++++++-
hysop/operator/spatial_filtering.py | 22 +--
9 files changed, 238 insertions(+), 56 deletions(-)
create mode 100644 hysop/backend/host/python/operator/spatial_filtering.py
diff --git a/examples/multiresolution/scalar_advection.py b/examples/multiresolution/scalar_advection.py
index db17fa924..7fc6e61c8 100644
--- a/examples/multiresolution/scalar_advection.py
+++ b/examples/multiresolution/scalar_advection.py
@@ -8,7 +8,7 @@ def compute(args):
from hysop.constants import Implementation
from hysop.operators import Advection, DirectionalAdvection, StrangSplitting, \
LowpassFilter, HDF_Writer
- from hysop.methods import Remesh, TimeIntegrator, ComputeGranularity, \
+ from hysop.methods import Remesh, TimeIntegrator, \
Interpolation
## Function to compute initial velocity values
@@ -22,19 +22,21 @@ def compute(args):
data[0][...] = args.velocity[::-1][i]
## Function to compute initial scalar values
- def init_scalar(data, coords):
- data[0][...] = 1.0
- for x in coords[0]:
- data[0][...] *= np.cos(x)
+ def init_scalar(data, coords, component=None):
+ if (component in (None,0)):
+ data[0][...] = 1.0
+ for x in coords[0]:
+ data[0][...] *= np.cos(x)
- data[1][...] = 1.0
- for x in coords[0]:
- data[1][...] *= np.sin(x)
+ if (component in (None,1)):
+ data[1][...] = 1.0
+ for x in coords[0]:
+ data[1][...] *= np.sin(x)
# Define domain
dim = args.ndim
npts = args.npts
- snpts = tuple(4*(x-1)+1 for x in args.npts)
+ snpts = tuple(4*x for x in args.npts)
box = Box(origin=args.box_origin, length=args.box_length, dim=dim)
# Get default MPI Parameters from domain (even for serial jobs)
@@ -49,7 +51,6 @@ def compute(args):
# Setup operator method dictionnary
# Advection-Remesh operator discretization parameters
method = {
- ComputeGranularity: args.compute_granularity,
TimeIntegrator: args.time_integrator,
Remesh: args.remesh_kernel,
Interpolation: args.interpolation
@@ -113,21 +114,27 @@ def compute(args):
# Finally insert our splitted advection into the problem
problem.insert(splitting)
+
+ #> Lowpass filters
+ lpfilter = LowpassFilter(input_variables={scalar: snpts},
+ output_variables={scalar: npts},
+ implementation=impl,
+ method={Remesh: Remesh.L2_1})
#> Operators to dump all fields
io_params = IOParams(filename='fine', frequency=args.dump_freq)
- df0 = HDF_Writer(name='S',
+ df0 = HDF_Writer(name='S_fin',
io_params=io_params,
variables={scalar: snpts},
**extra_op_kwds)
io_params = IOParams(filename='coarse', frequency=args.dump_freq)
- df1 = HDF_Writer(name='V',
+ df1 = HDF_Writer(name='S_coarse',
io_params=io_params,
- variables={velo: npts},
+ variables={scalar: npts},
**extra_op_kwds)
# Add a writer of input field at given frequency.
- problem.insert(df0, df1)
+ problem.insert(lpfilter, df0, df1)
problem.build(args)
# If a visu_rank was provided, and show_graph was set,
@@ -200,7 +207,7 @@ if __name__=='__main__':
parser = ScalarAdvectionArgParser()
parser.set_defaults(box_origin=(0.0,), box_length=(2*np.pi,),
- tstart=0.0, tend=2*np.pi, npts=(16,),
+ tstart=0.0, tend=2*np.pi, npts=(32,),
dump_freq=1, cfl=0.5, velocity=(1.0,),
ndim=3, compute_precision='fp64')
diff --git a/hysop/backend/host/python/operator/spatial_filtering.py b/hysop/backend/host/python/operator/spatial_filtering.py
new file mode 100644
index 000000000..13db58a20
--- /dev/null
+++ b/hysop/backend/host/python/operator/spatial_filtering.py
@@ -0,0 +1,27 @@
+
+from hysop.tools.types import check_instance, first_not_None
+from hysop.tools.decorators import debug
+from hysop.backend.host.host_operator import HostOperator
+from hysop.core.graph.graph import op_apply
+from hysop.fields.continuous_field import Field
+from hysop.parameters.parameter import Parameter
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.operator.base.spatial_filtering import LowpassFilterBase
+
+class PythonLowpassFilter(LowpassFilterBase, HostOperator):
+ """
+ Python implementation for lowpass spatial filtering: small grid -> coarse grid
+ """
+
+ @op_apply
+ def apply(self, **kwds):
+ """Apply analytic formula."""
+ super(PythonLowpassFilter, self).apply(**kwds)
+ fin, fout = self.fin, self.fout
+ iratio = self.iratio
+
+ fout[...] = 0
+ for (idx, Wi) in self.nz_weights.iteritems():
+ slc = tuple(slice(i, i+r*s, r) for (i,s,r) in zip(idx, fout.shape, iratio))
+ fout[...] += Wi*fin[slc]
+ self.dFout.exchange_ghosts()
diff --git a/hysop/core/graph/computational_node.py b/hysop/core/graph/computational_node.py
index 833cbe66e..adffa2b54 100644
--- a/hysop/core/graph/computational_node.py
+++ b/hysop/core/graph/computational_node.py
@@ -277,7 +277,7 @@ class ComputationalGraphNode(OperatorBase):
is_domainless = property(_get_is_domainless)
@classmethod
- def expand_tensor_fields(fields):
+ def expand_tensor_fields(cls, fields):
scalar_fields = ()
tensor_fields = ()
for field in fields:
diff --git a/hysop/core/graph/computational_operator.py b/hysop/core/graph/computational_operator.py
index 75dc455c2..f9953e7fb 100644
--- a/hysop/core/graph/computational_operator.py
+++ b/hysop/core/graph/computational_operator.py
@@ -58,6 +58,9 @@ class ComputationalGraphOperator(ComputationalGraphNode):
*override* the following method:
get_operator_requirements()
+ To declare that an operator does not invalidate an input_field when the
+ same field is written, *override* the get_preserved_input_fields method.
+
An operator has to *extend* the following abstract method:
apply(simulation, **kwds)
@@ -84,7 +87,7 @@ class ComputationalGraphOperator(ComputationalGraphNode):
-> topologies/topology descriptors are checked
initialize()
pre_initialize()
- handle method() -> self.method is set
+ handle method() -> self.method is set
get_field_requirements() -> self.input_field_requirements it set
get_operator_requirements() -> self.operator_requirements is set
check_requirements_compatibility() -> updates
@@ -98,7 +101,7 @@ class ComputationalGraphOperator(ComputationalGraphNode):
setup() (requires self.discretized to be set)
self.ready is set
apply() (requires self.ready to be set)
- abtract method to be implemented
+ abtract method to be implemented by each operator
finalize() (requires self.ready to be set)
sets self.ready to False
@@ -228,6 +231,14 @@ class ComputationalGraphOperator(ComputationalGraphNode):
requirements.update_outputs(output_field_requirements)
return requirements
+ def get_preserved_input_fields(self):
+ """
+ Declare fields that should not be invalidated as a set.
+ You can only declare scalar fields that are input and output fields at the same time.
+ For example, a redistribute operation does not invalidate input fields.
+ """
+ return set()
+
@debug
def get_node_requirements(self):
"""
@@ -427,6 +438,10 @@ class ComputationalGraphOperator(ComputationalGraphNode):
self.discrete_fields will be a tuple containing all input and
output discrete fields.
+
+ Discrete tensor fields are built back from discretized scalar fields
+ and are accessible from self.input_tensor_fields, self.output_tensor_fields
+ and self.discrete_tensor_fields like their scalar counterpart.
"""
if self.discretized:
return
diff --git a/hysop/core/graph/graph_builder.py b/hysop/core/graph/graph_builder.py
index 15365a3f6..d9930489b 100644
--- a/hysop/core/graph/graph_builder.py
+++ b/hysop/core/graph/graph_builder.py
@@ -369,8 +369,9 @@ class GraphBuilder(object):
istates = None if (current_level!=0) else input_states
cstate = self.topology_states.setdefault(ofield,
self.new_topology_state(ofield))
+ invalidate_field = (ofield not in op.get_preserved_input_fields())
dstate = cstate.handle_output(opnode, otopo, ofreqs,
- op, istates,
+ op, istates, invalidate_field,
graph, edge_properties, vertex_properties)
output_fields[ofield] = otopo
output_states[ofield] = dstate
@@ -711,6 +712,7 @@ class GraphBuilder(object):
# dictionnary (topology -> list of node) that are up to date (lastly written)
# multiple fields can be up to date at the same time after a redistribute operator
+ # or after an operator that implements the get_preserved_input_fields method.
self.write_nodes = {}
# dictionnary (topology -> list of nodes) that are currently reading field:topo
@@ -1096,7 +1098,7 @@ class GraphBuilder(object):
def handle_output(self, opnode, output_topo, oreqs,
- operator, input_topology_states,
+ operator, input_topology_states, invalidate_field,
graph, edge_properties, vertex_properties):
ofield = self.field
@@ -1112,19 +1114,19 @@ class GraphBuilder(object):
self.add_edge(graph, edge_properties,
write_nodes[output_topo], opnode,
ofield, output_topo)
-
- # add dependency to all operators that reads this field
- # to prevent concurent read-writes.
- if output_topo in read_nodes:
- for ro_node in read_nodes[output_topo]:
- self.add_edge(graph, edge_properties,
- ro_node, opnode,
- ofield, output_topo)
-
- # remove read/write dependencies and states
- # read_nodes.clear()
- write_nodes.clear()
- dtopology_states.clear()
+
+ if invalidate_field:
+ # add dependency to all operators that reads this field
+ # to prevent concurent read-writes.
+ if output_topo in read_nodes:
+ for ro_node in read_nodes[output_topo]:
+ self.add_edge(graph, edge_properties,
+ ro_node, opnode,
+ ofield, output_topo)
+
+ # remove read/write dependencies and states
+ write_nodes.clear()
+ dtopology_states.clear()
# add the operator node as the one that lastly wrote this field.
# no other operators can be reading as this topology just been written.
diff --git a/hysop/numerics/remesh/kernel_generator.py b/hysop/numerics/remesh/kernel_generator.py
index c7d1802de..364b83a32 100644
--- a/hysop/numerics/remesh/kernel_generator.py
+++ b/hysop/numerics/remesh/kernel_generator.py
@@ -150,16 +150,24 @@ class Kernel(object):
if verbose:
print ' => Generating lambdas'
-
+
+ imin=-Ms
+ imax=+Ms
if split_polys:
- gamma_l = sp.interpolate.PPoly(Cl,X,extrapolate=True)
- gamma_r = sp.interpolate.PPoly(Cr,X,extrapolate=True)
+ gamma_l = sp.interpolate.PPoly(Cl,X,extrapolate=False)
+ gamma_r = sp.interpolate.PPoly(Cr,X,extrapolate=False)
def gamma(x):
i = np.floor(x)
z = x-i
- return (z>=0.5)*gamma_r(i+1-z) + (z<0.5)*gamma_l(x)
+ res = (z>=0.5)*gamma_r(i+1-z) + (z<0.5)*gamma_l(x)
+ res[np.isnan(res)] = 0.0
+ return res
else:
- gamma = sp.interpolate.PPoly(C,X,extrapolate=True)
+ gamma_ = sp.interpolate.PPoly(C,X,extrapolate=False)
+ def gamma(x):
+ res = gamma_(x)
+ res[np.isnan(res)] = 0.0
+ return res
self.I = X
self.gamma = gamma
@@ -179,7 +187,9 @@ class Kernel(object):
def __call__(self,x):
Ms = self.Ms
- return self.gamma(x)
+ res = self.gamma(x)
+ return res
+
class SymmetricKernelGenerator(object):
@@ -433,7 +443,7 @@ if __name__=='__main__':
fig = plt.figure()
plt.xlabel(r'$x$')
plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p,'r')+'}$')
- X = np.linspace(-k0.Ms,+k0.Ms,1000)
+ X = np.linspace(-k0.Ms-1,+k0.Ms+1,1000)
s = plt.subplot(1,1,1)
for i,k in enumerate(kernels):
s.plot(X,k(X),label=r'$\Lambda_{'+'{},{}'.format(p,k.r)+'}$')
diff --git a/hysop/numerics/remesh/remesh.py b/hysop/numerics/remesh/remesh.py
index ad7da6671..bf2c38bd1 100644
--- a/hysop/numerics/remesh/remesh.py
+++ b/hysop/numerics/remesh/remesh.py
@@ -73,7 +73,7 @@ if __name__=='__main__':
fig = plt.figure()
plt.xlabel(r'$x$')
plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p,'r')+'}$')
- X = np.linspace(-k0.Ms,+k0.Ms,1000)
+ X = np.linspace(-k0.Ms-1,+k0.Ms+1,1000)
s = plt.subplot(1,1,1)
for i,k in enumerate(kernels):
s.plot(X,k(X),label=r'$\Lambda_{'+'{},{}'.format(p,k.r)+'}$')
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index f160158e5..859b68a52 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -3,13 +3,16 @@
LowpassFilter operator generator.
"""
from hysop.constants import Implementation
-from hysop.tools.types import check_instance, to_list, first_not_None
+from hysop.tools.types import check_instance, to_list, first_not_None, InstanceOf
+from hysop.tools.numpywrappers import npw
from hysop.tools.decorators import debug
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import Field, ScalarField
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.core.graph.node_generator import ComputationalGraphNodeGenerator
-from hysop.core.graph.computational_node_fronted import ComputationalGraphNodeFrontend
+from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
+from hysop.methods import Remesh
+from hysop.numerics.remesh.remesh import RemeshKernel
class LowpassFilterBase(object):
@@ -17,6 +20,15 @@ class LowpassFilterBase(object):
Base implementation for lowpass spatial filtering: small grid -> coarse grid
"""
+ __default_method = {
+ Remesh: Remesh.L2_1,
+ }
+
+ __available_methods = {
+ Remesh: (InstanceOf(Remesh), InstanceOf(RemeshKernel)),
+ }
+
+
def __init__(self, input_field, output_field,
input_topo, output_topo,
**kwds):
@@ -33,7 +45,108 @@ class LowpassFilterBase(object):
self.Fin = input_field
self.Fout = output_field
+
+ @classmethod
+ def default_method(cls):
+ dm = super(LowpassFilterBase, cls).default_method()
+ dm.update(cls.__default_method)
+ return dm
+
+ @classmethod
+ def available_methods(cls):
+ am = super(LowpassFilterBase, cls).available_methods()
+ am.update(cls.__available_methods)
+ return am
+
+ @classmethod
+ def supports_multiple_field_topologies(cls):
+ return True
+
+ @classmethod
+ def supports_mpi(cls):
+ return True
+
+ def get_preserved_input_fields(self):
+ return {self.Fout}
+
+ @debug
+ def handle_method(self,method):
+ super(LowpassFilterBase, self).handle_method(method)
+ remesh_kernel = method.pop(Remesh)
+ if isinstance(remesh_kernel, Remesh):
+ remesh_kernel = RemeshKernel.from_enum(remesh_kernel)
+ self.remesh_kernel = remesh_kernel
+
+ @classmethod
+ def remesh_ghosts(cls, remesh_kernel):
+ """Return the minimum number of ghosts for remeshed scalars."""
+ assert remesh_kernel.n >=1 , 'Bad remeshing kernel.'
+ if remesh_kernel.n >1:
+ assert remesh_kernel.n % 2 == 0, 'Odd remeshing kernel moments.'
+ min_ghosts = int(remesh_kernel.n//2)+1
+ return min_ghosts
+
+ @debug
+ def get_field_requirements(self):
+ requirements = super(LowpassFilterBase, self).get_field_requirements()
+ dim = self.Fin.dim
+
+ Fin_topo, Fin_requirements = requirements.get_input_requirement(self.Fin)
+ try:
+ Fin_dx = Fin_topo.space_step
+ except AttributeError:
+ Fin_dx = Fin_topo.mesh.space_step
+
+ Fout_topo, Fout_requirements = requirements.get_output_requirement(self.Fout)
+ try:
+ Fout_dx = Fout_topo.space_step
+ except AttributeError:
+ Fout_dx = Fout_topo.mesh.space_step
+
+ ratio = Fout_dx / Fin_dx
+ msg='Destination grid is finer than source grid: {}'.format(ratio)
+ assert (ratio>=1.0).all(), msg
+ iratio = ratio.astype(npw.int32)
+ msg='Grid ratio is not an integer on at least one axis: {}'.format(ratio)
+ assert (ratio==iratio).all(), msg
+
+ remesh_ghosts = self.remesh_ghosts(self.remesh_kernel)
+ fine_grid_ghosts = iratio*remesh_ghosts - 1
+ Fin_requirements.min_ghosts = fine_grid_ghosts
+
+ self.remesh_ghosts = remesh_ghosts
+ self.fine_grid_ghosts = fine_grid_ghosts
+
+ return requirements
+
+ def compute_weights(self, iratio, product=True):
+ assert (iratio>=1).all()
+ remesh_kernel = self.remesh_kernel
+ p = remesh_kernel.n//2 + 1
+ shape = 2*p*iratio-1
+ weights = npw.zeros(dtype=npw.float64, shape=shape)
+ nz_weights = {}
+ for idx in npw.ndindex(*shape):
+ X = (npw.asarray(idx, dtype=npw.float64)+1) / iratio - p
+ if product:
+ W = npw.prod(remesh_kernel(X))
+ else:
+ # this does not seem to work
+ R = npw.sqrt(npw.dot(X,X))
+ W = remesh_kernel(R)
+ weights[idx] = W
+ if (W!=0):
+ nz_weights[idx] = W
+ Ws = weights.sum()
+ weights = weights / Ws
+ nz_weights = {k: v/Ws for (k,v) in nz_weights.iteritems()}
+
+ assert abs(weights.sum() - 1.0) < 1e-8, weights.sum()
+ assert abs(npw.sum(nz_weights.values()) - 1.0) < 1e-8, npw.sum(nz_weights.values())
+ self.iratio = iratio
+ self.weights = weights
+ self.nz_weights = nz_weights
@debug
def discretize(self):
@@ -43,10 +156,16 @@ class LowpassFilterBase(object):
dFin = self.get_input_discrete_field(self.Fin)
dFout = self.get_output_discrete_field(self.Fout)
- import sys
- sys.exit(1)
+ iratio = dFin.compute_resolution / dFout.compute_resolution
+ self.compute_weights(iratio)
+
+ remesh_ghosts = self.remesh_ghosts
+ fine_grid_ghosts = iratio*remesh_ghosts - 1
+ fin = dFin.sdata[dFin.local_slices(ghosts=fine_grid_ghosts)]
+ fout = dFout.compute_buffers[0]
self.dFin, self.dFout = dFin, dFout
+ self.fin, self.fout = fin, fout
diff --git a/hysop/operator/spatial_filtering.py b/hysop/operator/spatial_filtering.py
index 7572dc482..feccbf392 100644
--- a/hysop/operator/spatial_filtering.py
+++ b/hysop/operator/spatial_filtering.py
@@ -5,9 +5,10 @@ LowpassFilter operator generator.
from hysop.constants import Implementation
from hysop.tools.types import check_instance, to_list, first_not_None
from hysop.tools.decorators import debug
-from hysop.fields.continuous_field import Field
+from hysop.fields.continuous_field import Field, ScalarField
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.core.graph.computational_node import ComputationalGraphNode
from hysop.core.graph.node_generator import ComputationalGraphNodeGenerator
from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
@@ -66,7 +67,7 @@ class LowpassFilterFrontend(ComputationalGraphNodeFrontend):
check_instance(output_field, ScalarField)
check_instance(input_topo, CartesianTopologyDescriptors)
check_instance(output_topo, CartesianTopologyDescriptors)
- check_instance(base_kwds, dict, keys=str)
+ check_instance(base_kwds, dict, keys=str, allow_none=True)
assert (input_topo != output_topo), "Same topology for input and output."
super(LowpassFilterFrontend, self).__init__(input_field=input_field, input_topo=input_topo,
@@ -80,7 +81,7 @@ class LowpassFilter(ComputationalGraphNodeGenerator):
Graphnode generator to lowpass filter multiple fields at once.
"""
@debug
- def __init__(self, input_variables, output_variables=None,
+ def __init__(self, input_variables, output_variables,
implementation=None,
base_kwds=None,
**kwds):
@@ -104,11 +105,9 @@ class LowpassFilter(ComputationalGraphNodeGenerator):
kwds:
Extra parameters passed to generated operators.
"""
- input_fields = list(self.expand_tensor_fields(input_variables.keys())[0])
+ input_fields = list(ComputationalGraphNode.expand_tensor_fields(input_variables.keys())[0])
assert len(set(input_fields)) == len(input_fields)
- if (output_variables is None):
- output_variables = input_variables
- output_fields = list(self.expand_tensor_fields(output_variables.keys())[0])
+ output_fields = list(ComputationalGraphNode.expand_tensor_fields(output_variables.keys())[0])
assert len(input_fields) == len(output_fields)
for i, field in enumerate(output_fields):
if (field is None):
@@ -124,7 +123,10 @@ class LowpassFilter(ComputationalGraphNodeGenerator):
check_instance(output_variables, dict, keys=Field, values=CartesianTopologyDescriptors)
check_instance(base_kwds, dict, keys=str)
- super(LowpassFilter, self).__init__(**base_kwds)
+ super(LowpassFilter, self).__init__(
+ candidate_input_tensors=None,
+ candidate_output_tensors=None,
+ **base_kwds)
self._input_fields = input_fields
self._output_fields = output_fields
@@ -137,8 +139,8 @@ class LowpassFilter(ComputationalGraphNodeGenerator):
def _generate(self):
nodes = []
for (ifield, ofield) in zip(self._input_fields, self._output_fields):
- stopo = self.get_topo_descriptor(self._input_variables, ifield)
- ttopo = self.get_topo_descriptor(self._output_variables, ofield)
+ stopo = ComputationalGraphNode.get_topo_descriptor(self._input_variables, ifield)
+ ttopo = ComputationalGraphNode.get_topo_descriptor(self._output_variables, ofield)
impl = self._impl
kwds = self._kwds.copy()
--
GitLab