diff --git a/examples/cpu_plane_jet/cpu_planejet.py b/examples/cpu_plane_jet/cpu_planejet.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ba17a6a7034f75e3dbc53642b093aa11dd87daf
--- /dev/null
+++ b/examples/cpu_plane_jet/cpu_planejet.py
@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+
+import hysop
+
+from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh,\
+ Support, Splitting, MultiScale, MultiScale, SpaceDiscretisation, \
+ GhostUpdate, Scales, dtCrit, ExtraArgs
+
+from hysop.problem.simulation import Simulation
+from hysop.fields.variable_parameter import VariableParameter
+
+from hysop.numerics.odesolvers import RK2 as RK2
+from hysop.numerics.odesolvers import RK3 as RK3
+from hysop.numerics.odesolvers import RK4 as RK4
+
+from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.remeshing import L6_4 as remesh_formula
+from hysop.numerics.remeshing import Linear
+
+from hysop.operator.advection import Advection
+from hysop.operator.diffusion import Diffusion
+from hysop.operator.stretching import Stretching
+from hysop.operator.poisson import Poisson
+from hysop.operator.adapt_timestep import AdaptTimeStep
+
+from hysop.operator.redistribute_inter import RedistributeInter
+from hysop.operator.redistribute_intra import RedistributeIntra
+
+from hysop.gpu.gpu_transfer import DataTransfer
+from hysop.operator.hdf_io import HDF_Writer
+
+import utilities as utils
+
+
+## Simulation parameters
+dim = 3
+N = 1<<8 + 1
+nb_elems = [N]*dim
+
+## Physical parameters:
+# Flow viscosity
+viscosity = 1e-4
+
+## I/O parameters
+dt_output = IOParams(frequency=out_freq, filename='cpu_dt.dat', fileformat=ASCII)
+velo_output = IOParams(frequency=out_freq, filename='cpu_velo.dat', fileformat=HDF5)
+vorti_output = IOParams(frequency=out_freq, filename='cpu_vorti.dat', fileformat=HDF5)
+
+
+# Domain
+box = hysop.Box(length=[1.0]*dim,
+ origin=[0.0]*dim,
+ bc=PERIODIC)
+
+# Fields
+velo = hysop.Field(domain=box, formula=utils.initialize_velocity,
+ name='Velocity', is_vector=True)
+vorti = hysop.Field(domain=box,
+ name='Vorticity', is_vector=True)
+
+# Discretizations
+d_uw = Discretization(nb_elems)
+
+# Topologies
+topo_uw = box.create_topology(d_uw, dim=1)
+
+# Adaptative timestep simulation
+dt = VariableParameter({'dt': 0.001})
+simu = Simulation(start=0.0, end=5.0, time_step=0.001, max_iter=10000)
+
+
+## Operators
+ops = {}
+
+# CPU operators
+ops['advec_w'] = Advection(velo,
+ discretization=topo_uw,
+ variables={vorti: topo_uw},
+ method={Scales: 'p_64', MultiScale: 'L4_4'})
+ops['stretching'] = Stretching(velo, vorti, discretization=d_uw)
+ops['diffusion'] = Diffusion(variables={vorti: d_uw}, viscosity=viscosity)
+ops['poisson'] = Poisson(velo, vorti, discretization=d_uw)
+
+
+# other operators
+ops['dt_adapt'] = AdaptTimeStep(velo, vorti,
+ simulation=simu,
+ time_range=[0, np.infty],
+ discretization=d_uw,
+ method={TimeIntegrator: RK3,
+ SpaceDiscretisation: FD_C_4,
+ dtCrit: ['gradU', 'cfl']},
+ lcfl=0.15,
+ cfl=1.5,
+ io_params=dt_output)
+
+ops['velo_io'] = HDF_Writer(variables={velo : topo_uw}, io_params=velo_output)
+ops['vorti_io'] = HDF_Writer(variables={vorti: topo_uw}, io_params=vorti_output)
+
+for op in ops:
+ op.discretize()
+
+
diff --git a/examples/cpu_plane_jet/create_random_arrays.py b/examples/cpu_plane_jet/create_random_arrays.py
new file mode 100644
index 0000000000000000000000000000000000000000..de7322631e99cba53ce5e92b837ad3074bf21349
--- /dev/null
+++ b/examples/cpu_plane_jet/create_random_arrays.py
@@ -0,0 +1,39 @@
+import os
+import numpy as np
+from hysop.constants import HYSOP_REAL, ORDER
+
+
+def random_init(shape, mpi_comm):
+ # Create a folder to store all random arrays
+ d = 'rand_init'
+ if mpi_comm.Get_rank() == 0:
+ if not os.path.exists(d):
+ os.makedirs(d)
+ mpi_comm.Barrier()
+ file_name = "{0}_{1}_{2}".format(*shape)
+ file_name += "_{0}p_{1}.dat".format(mpi_comm.Get_size(),
+ mpi_comm.Get_rank())
+ try:
+ randX = np.asarray(
+ np.reshape(np.fromfile(os.path.join(d, 'randX_' + file_name)),
+ shape),
+ dtype=HYSOP_REAL, order=ORDER)
+ randY = np.asarray(
+ np.reshape(np.fromfile(os.path.join(d, 'randY_' + file_name)),
+ shape),
+ dtype=HYSOP_REAL, order=ORDER)
+ randZ = np.asarray(
+ np.reshape(np.fromfile(os.path.join(d, 'randZ_' + file_name)),
+ shape),
+ dtype=HYSOP_REAL, order=ORDER)
+ except IOError:
+ randX = np.asarray(np.random.random(shape),
+ dtype=HYSOP_REAL, order=ORDER) - 0.5
+ randY = np.asarray(np.random.random(shape),
+ dtype=HYSOP_REAL, order=ORDER) - 0.5
+ randZ = np.asarray(np.random.random(shape),
+ dtype=HYSOP_REAL, order=ORDER) - 0.5
+ randX.tofile(os.path.join(d, 'randX_' + file_name))
+ randY.tofile(os.path.join(d, 'randY_' + file_name))
+ randZ.tofile(os.path.join(d, 'randZ_' + file_name))
+ return randX, randY, randZ
diff --git a/examples/cpu_plane_jet/utilities.py b/examples/cpu_plane_jet/utilities.py
new file mode 100644
index 0000000000000000000000000000000000000000..f09774fe55ee70b011436fc66e40577a4b1e474c
--- /dev/null
+++ b/examples/cpu_plane_jet/utilities.py
@@ -0,0 +1,27 @@
+
+import numpy as np
+
+width = 0.01
+ampl1 = 0.05
+ampl3 = 0.3
+
+
+def initialize_velocity(res, x, y, z, t, rand_init=False):
+ yy = np.abs(y - 0.5)
+ aux = (0.1 - 2.0 * yy) / (4.0 * width)
+ strg = np.exp(np.abs(aux**2))
+ if rand_init:
+ from create_random_arrays import random_init
+ randX, randY, randZ = random_init(res[0].shape, box.comm_task)
+ strg1 = strg * randX
+ strg2 = strg * randY
+ strg3 = strg * randZ
+ res[0][...] = 0.5 * (1.0 + np.tanh(aux)) * (1.0 + ampl3 * np.sin(8.0 * np.pi * x)) * (1.0 + ampl1 * strg1)
+ res[1][...] = ampl1 * strg2
+ res[2][...] = ampl1 * strg3
+ else:
+ res[0][...] = 0.5 * (1.0 + np.tanh(aux)) * (1.0 + ampl3 * np.sin(8.0 * np.pi * x))
+ res[1][...] = 0.0
+ res[2][...] = 0.0
+ return res
+