diff --git a/examples/example_utils.py b/examples/example_utils.py
index 46cbc096c1ad778eb0201d38e10e52fc41bab144..879fed10716cea515bc925858f23ddec7a6f7234 100644
--- a/examples/example_utils.py
+++ b/examples/example_utils.py
@@ -553,6 +553,10 @@ class HysopArgParser(argparse.ArgumentParser):
default='LINEAR',
dest='polynomial_interpolator',
help='Set the default polynomial interpolator for polynomial interpolation or restriction methods.')
+ method.add_argument('-ai', '--advection-interpolator', type=str,
+ default='LEGACY',
+ dest='advection_interpolator',
+ help='Set the default polynomial interpolator for bilevel advection. Legacy interpolator uses linear interpolation, else specify a custom polynomial interpolator.')
method.add_argument('-sf', '--stretching-formulation', type=str,
default='conservative',
dest='stretching_formulation',
@@ -594,6 +598,8 @@ class HysopArgParser(argparse.ArgumentParser):
allow_subgrid=False)
args.polynomial_interpolator = \
self._convert_polynomial_interpolation('polynomial_interpolator', args.polynomial_interpolator)
+ args.advection_interpolator = \
+ self._convert_advection_interpolation('advection_interpolator', args.advection_interpolator)
def _add_threading_args(self):
threading = self.add_argument_group('threading parameters')
@@ -1355,6 +1361,27 @@ class HysopArgParser(argparse.ArgumentParser):
}
return self._check_convert(argname, pi, polynomial_interpolations)
+ def _convert_advection_interpolation(self, argname, ai):
+ from hysop.numerics.interpolation.interpolation import Interpolation
+ from hysop.numerics.interpolation.polynomial import PolynomialInterpolation
+ advection_interpolations = {
+ 'legacy': Interpolation.LINEAR,
+ 'linear': PolynomialInterpolation.LINEAR,
+ 'cubic_fdc2': PolynomialInterpolation.CUBIC_FDC2,
+ 'cubic_fdc4': PolynomialInterpolation.CUBIC_FDC4,
+ 'cubic_fdc6': PolynomialInterpolation.CUBIC_FDC6,
+ 'quintic_fdc2': PolynomialInterpolation.QUINTIC_FDC2,
+ 'quintic_fdc4': PolynomialInterpolation.QUINTIC_FDC4,
+ 'quintic_fdc6': PolynomialInterpolation.QUINTIC_FDC6,
+ 'septic_fdc2': PolynomialInterpolation.SEPTIC_FDC2,
+ 'septic_fdc4': PolynomialInterpolation.SEPTIC_FDC4,
+ 'septic_fdc6': PolynomialInterpolation.SEPTIC_FDC6,
+ 'nonic_fdc2': PolynomialInterpolation.NONIC_FDC2,
+ 'nonic_fdc4': PolynomialInterpolation.NONIC_FDC4,
+ 'nonic_fdc6': PolynomialInterpolation.NONIC_FDC6,
+ }
+ return self._check_convert(argname, ai, advection_interpolations)
+
def _convert_stretching_formulation(self, argname, stretching_formulation):
from hysop.constants import StretchingFormulation
stretching_formulations = {
diff --git a/examples/scalar_advection/levelset.py b/examples/scalar_advection/levelset.py
index c455ef1be2b0e9f67a7562452114bf082c012e36..28797455179160f5ba50f2d81393161efa90da60 100644
--- a/examples/scalar_advection/levelset.py
+++ b/examples/scalar_advection/levelset.py
@@ -73,8 +73,9 @@ def compute(args):
# Setup operator method dictionnary
method.update({
- TimeIntegrator: args.time_integrator,
- Remesh: args.remesh_kernel,
+ TimeIntegrator: args.time_integrator,
+ Remesh: args.remesh_kernel,
+ Interpolation: args.advection_interpolator
})
# Create a simulation and solve the problem
diff --git a/hysop/backend/device/codegen/kernels/directional_advection.py b/hysop/backend/device/codegen/kernels/directional_advection.py
index aba77d14c2193dd8df3bd01977073706836c78a1..576be34f6d5bc838fec9eb29b13709645e480d9d 100644
--- a/hysop/backend/device/codegen/kernels/directional_advection.py
+++ b/hysop/backend/device/codegen/kernels/directional_advection.py
@@ -102,7 +102,6 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
if (polynomial_interpolator is not None):
if is_bilevel:
- print 'LOL'
assert (integer_grid_ratio == grid_ratio).all()
assert (integer_grid_ratio[::-1] == polynomial_interpolator.gr).all()
else:
@@ -185,8 +184,6 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
self.polynomial_interpolator = polynomial_interpolator
self.gencode()
- self.edit()
-
@classmethod
def advection_ghosts(cls, velocity_cfl):
@@ -202,11 +199,13 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
return linear_interpolation_ghosts+int(math.floor(grid_ratio*velocity_cfl))
@classmethod
- def min_velocity_ghosts(cls, work_dim, velocity_cfl, grid_ratio, periodicity):
+ def min_velocity_ghosts(cls, work_dim, velocity_cfl, grid_ratio, periodicity, polynomial_interpolator):
"""Return the minimal numbers of ghosts required on the velocity grid, includint interpolation ghosts."""
interpolation_ghosts = 1
ghosts = periodicity * (grid_ratio > 1.0)
ghosts += interpolation_ghosts
+ if (polynomial_interpolator is not None):
+ ghosts += polynomial_interpolator.ghosts
ghosts[-1] += cls.advection_ghosts(velocity_cfl)
return np.asarray(ghosts, dtype=np.int32)
@@ -491,7 +490,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
s.append('{} = {} - {};'.format(
velocity_h[i], velocity_p[i], velocity_gid[i]))
if polynomial_interpolator:
- s.append('{} = (int)(round({}/{}));'.format(widx[i], velocity_h[i], p_dx[i]))
+ s.append('{} = (int)(round({}*{}));'.format(widx[i], velocity_h[i], grid_ratio[i]))
else:
s.append('{} = {} + {};'.format(velocity_gid[i], 'kji'[i], velocity_ghosts))
yield ctx
@@ -585,7 +584,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
velocity_h[i], velocity_p[i], velocity_gid[i]))
is_first = True
if polynomial_interpolator:
- s.append('{} = (int)(round({}/{}));'.format(widx[i], velocity_h[i], p_dx[i]))
+ s.append('{} = (int)(round({}*{}));'.format(widx[i], velocity_h[i], grid_ratio[i]))
for idx in np.ndindex(*polynomial_interpolator.n):
idx = idx[::-1]
offset = '+'.join('({}{:+})*{}'.format(gidi, idxi-gi, vsi)
diff --git a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
index f5268ea2c88a551df2bc0a71c6a7e7ab18b081b4..a5b19a0da0ed71988ac994b26358a2a7ff7add69 100644
--- a/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
+++ b/hysop/backend/device/opencl/autotunable_kernels/advection_dir.py
@@ -33,14 +33,13 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
raise ValueError(msg)
if not isinstance(velocity_cfl, float) or velocity_cfl<=0.0:
raise ValueError('Invalid velocity_cfl value.')
- check_instance(polynomial_interpolator, PolynomialSubgridInterpolator)
+ check_instance(polynomial_interpolator, PolynomialSubgridInterpolator, allow_none=True)
check_instance(grid_ratio, npw.ndarray, minval=1.0, allow_none=True)
is_bilevel = (grid_ratio>1.0).any()
periodicity = velocity.periodicity.astype(npw.int32)
-
- min_velocity_ghosts = DirectionalAdvectionKernelGenerator.min_velocity_ghosts(dim, velocity_cfl, grid_ratio, periodicity)
+ min_velocity_ghosts = DirectionalAdvectionKernelGenerator.min_velocity_ghosts(dim, velocity_cfl, grid_ratio, periodicity, polynomial_interpolator)
velocity_cache_ghosts = DirectionalAdvectionKernelGenerator.advection_cache_ghosts(velocity_cfl, grid_ratio[-1])
self.check_cartesian_field(velocity, min_ghosts=min_velocity_ghosts, nb_components=1)
@@ -60,10 +59,19 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
Vr = relative_velocity
assert isinstance(Vr, (str,float))
+
+ if is_bilevel:
+ if (polynomial_interpolator is not None):
+ sbilevel = '_bilevel_{}'.format(
+ getattr(polynomial_interpolator.interpolator, 'spi', 'unknown_interpolator'))
+ else:
+ sbilevel='_bilevel_legacy_linear'
+ else:
+ sbilevel = ''
ftype = clTools.dtype_to_ctype(precision)
name = 'directional_advection{}_{}__{}_{}__{}'.format(
- "" if (is_bilevel is None) else "_bilevel_linear",
+ sbilevel,
DirectionLabels[direction],
time_integrator.name(), ftype, abs(hash(Vr)))
@@ -241,7 +249,8 @@ class OpenClAutotunableDirectionalAdvectionKernel(OpenClAutotunableKernel):
def hash_extra_kwds(self, extra_kwds):
"""Hash extra_kwds dictionnary for caching purposes."""
kwds = ('rk_scheme', 'ftype', 'work_dim', 'Vr', 'grid_ratio',
- 'vboundaries', 'cache_ghosts', 'work_size', 'known_args')
+ 'vboundaries', 'cache_ghosts', 'work_size', 'known_args',
+ 'polynomial_interpolator')
return self.custom_hash(*tuple(extra_kwds[kwd] for kwd in kwds),
mesh_info_vars=extra_kwds['mesh_info_vars'])
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index 2a2584eed5fbe48a30e0340f5c3563f23dc9622b..c5f969119cd48da468607e73a1812973de6e7b7b 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -187,6 +187,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
lidx: int32 empty input buffer, destroyed
ridx: int32 empty input buffer, destroyed
"""
+ assert self.interp is Interpolation.LINEAR # legacy linear interpolator
N = dX.shape[-1]
Iy = I[:-1]
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index 486c476285cce4f8aebaf3304d2562dd9596ca6e..d93bbefefb1d46b0ca01a3b22d2f378115b539b0 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -81,7 +81,9 @@ class PolynomialInterpolator(object):
if (fd is None):
msg='Could not determine finite differences order.'
raise RuntimeError(msg)
- return cls(deg=deg, fd=fd, **kwds)
+ obj = cls(deg=deg, fd=fd, **kwds)
+ obj.spi = spi
+ return obj
@classmethod
def cache_file(cls):
@@ -415,6 +417,10 @@ class PolynomialInterpolator(object):
def generate_subgrid_interpolator(self, grid_ratio, dtype=None):
return PolynomialSubgridInterpolator(interpolator=self,
grid_ratio=grid_ratio, dtype=dtype)
+
+ def __hash__(self):
+ objs = (self.dim, self.deg, self.fd, self.approximative)
+ return hash(objs)
class PolynomialSubgridInterpolator(object):
@@ -602,6 +608,10 @@ class PolynomialSubgridInterpolator(object):
Wr = self.W.reshape(self.s+self.n)[view].reshape(self.GR, self.N).copy()
return Wr
+ def __hash__(self):
+ objs = (self.interpolator, self.gr)
+ return hash(objs)
+
class PolynomialSubgridRestrictor(object):
def __init__(self, subgrid_interpolator):
diff --git a/hysop/operator/tests/test_bilevel_advection.py b/hysop/operator/tests/test_bilevel_advection.py
index d72b27dddba5bf5e4bf63dd4af2276bb8898b6e5..8eb2867d505cef14bb687f98a1a520864c02fe03 100644
--- a/hysop/operator/tests/test_bilevel_advection.py
+++ b/hysop/operator/tests/test_bilevel_advection.py
@@ -12,12 +12,13 @@ from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.operator.advection import Advection
from hysop import Field, Box
-from hysop.methods import Remesh, TimeIntegrator
+from hysop.methods import Remesh, TimeIntegrator, Interpolation
from hysop.constants import Implementation, DirectionLabels, Backend, \
HYSOP_REAL, Implementation
from hysop.numerics.splitting.strang import StrangSplitting, StrangOrder
from hysop.numerics.odesolvers.runge_kutta import Euler, RK2, RK4
from hysop.numerics.remesh.remesh import RemeshKernel
+from hysop.numerics.interpolation.polynomial import PolynomialInterpolation
class TestBilevelAdvectionOperator(object):
@@ -119,29 +120,37 @@ class TestBilevelAdvectionOperator(object):
implementations.remove(Implementation.PYTHON) # no bilevel support in python
- method = {TimeIntegrator: time_integrator, Remesh: remesh_kernel}
+ method = {
+ TimeIntegrator: time_integrator,
+ Remesh: remesh_kernel,
+ }
- def iter_impl(impl):
+ def iter_impl(impl, method=method):
graph = ComputationalGraph(name='test_graph')
+ method[Interpolation] = Interpolation.LINEAR
if (impl is Implementation.OPENCL):
- for cl_env in iter_clenv():
- msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
- cl_env.device.name.strip())
- da = DirectionalAdvection(
- velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,
- velocity_cfl=velocity_cfl, variables=variables, implementation=impl,
- method=method, name='advection_{}'.format(str(impl).lower()))
- split = StrangSplitting(
- splitting_dim=dim,
- extra_kwds=dict(cl_env=cl_env),
- order=StrangOrder.STRANG_SECOND_ORDER)
- force_tstate = ForceTopologyState(fields=variables.keys(),
- variables=variables,
- backend=Backend.OPENCL,
- extra_kwds={'cl_env': cl_env})
- split.push_operators(da)
- graph.push_nodes(split, force_tstate)
- yield msg, graph
+ for interp_method in (Interpolation.LINEAR, PolynomialInterpolation.LINEAR):
+ graph = ComputationalGraph(name='test_graph')
+ method[Interpolation] = interp_method
+ for cl_env in iter_clenv():
+ msg='platform {}, device {}, {}::{}'.format(cl_env.platform.name.strip(),
+ cl_env.device.name.strip(),
+ type(interp_method), interp_method)
+ da = DirectionalAdvection(velocity=vin,
+ advected_fields=sin, advected_fields_out=sout, dt=dt,
+ velocity_cfl=velocity_cfl, variables=variables, implementation=impl,
+ method=method, name='advection_{}'.format(str(impl).lower()))
+ split = StrangSplitting(
+ splitting_dim=dim,
+ extra_kwds=dict(cl_env=cl_env),
+ order=StrangOrder.STRANG_SECOND_ORDER)
+ force_tstate = ForceTopologyState(fields=variables.keys(),
+ variables=variables,
+ backend=Backend.OPENCL,
+ extra_kwds={'cl_env': cl_env})
+ split.push_operators(da)
+ graph.push_nodes(split, force_tstate)
+ yield msg, graph
elif impl is Implementation.FORTRAN:
assert dim==3, "Scales is only 3D"
adv = Advection(velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,