diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
deleted file mode 100755
index 54a4dadfb39bfbffdfbca2a00bbded0254c65375..0000000000000000000000000000000000000000
--- a/Examples/NSDebug.py
+++ /dev/null
@@ -1,472 +0,0 @@
-#!/usr/bin/python
-
-"""
-Navier Stokes 3D : flow past bluff bodies (using penalization).
-
-All parameters are set and defined in python module dataNS_bb.
-
-"""
-
-import parmepy as pp
-from parmepy.f2py import fftw2py
-import numpy as np
-from parmepy.problem.simulation import Simulation
-from parmepy.domain.obstacle.controlBox import ControlBox
-from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
- Remesh, Support, Splitting, dtCrit, SpaceDiscretisation, GhostUpdate
-from parmepy.numerics.integrators.runge_kutta2 import RK2 as RK2
-from parmepy.numerics.integrators.runge_kutta3 import RK3 as RK3
-from parmepy.numerics.integrators.runge_kutta4 import RK4 as RK4
-from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
-from parmepy.numerics.interpolation import Linear
-from parmepy.numerics.remeshing import L6_4 as rmsh
-import parmepy.tools.numpywrappers as npw
-from parmepy.constants import HDF5
-import parmepy.tools.io_utils as io
-
-
-## ----------- A 3d problem -----------
-print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
-
-## pi constant
-pi = np.pi
-cos = np.cos
-sin = np.sin
-
-# ====== Flow constants =======
-uinf = 1.0
-VISCOSITY = 1. / 1600.
-
-# ========= Geometry of the domain =========
-dim = 3
-#boxlength = npw.realarray([2.0, 2.0, 2.0])
-#boxorigin = npw.realarray([-1., -1., -1.])
-boxlength = npw.realarray([3., 5.12, 5.12])
-boxorigin = npw.realarray([-2.0, -2.56, -2.56])
-# The domain
-box = pp.Box(dim, length=boxlength, origin=boxorigin)
-
-# Sphere inside the domain
-RADIUS = 0.5
-sphere_pos = [0., 0., 0.]
-from parmepy.domain.obstacle.sphere import Sphere
-sphere = Sphere(domain=box, position=sphere_pos, radius=RADIUS)
-
-# ========= Discretisation parameters =========
-#nbElem = [1025, 513, 513]
-nbElem = [129, 65, 65]
-#nb = 17
-#nbElem = [nb, nb, nb]
-
-# ========= Vector Fields and variable parameters =========
-
-# Function to compute initial velocity
-def computeVel(res, x, y, z, t):
- ## res[0][...] = sin(x) * cos(y) * cos(z)
- ## res[1][...] = - cos(x) * sin(y) * cos(z)
- ## res[2][...] = 0.
- res[0][...] = uinf
- res[1][...] = 0.
- res[2][...] = 0.
- return res
-
-# Function to compute initial vorticity
-def computeVort(res, x, y, z, t):
- res[0][...] = 0.
- res[1][...] = 0.
- res[2][...] = 0.
- return res
-
-from parmepy import VariableParameter, Field
-velo = Field(domain=box, formula=computeVel,
- name='Velocity', isVector=True)
-vorti = Field(domain=box, formula=computeVort,
- name='Vorticity', isVector=True)
-dt = VariableParameter(data=0.0125, name='dt')
-
-# ========= Operators that describe the problem =========
-op = {}
-
-# Adaptative timestep method : dt = min(values(dtCrit))
-# where dtCrit is a list of criterions on which the computation
-# of the adaptative time step is based
-# ex : dtCrit = ['gradU', 'cfl', 'stretch'], means :
-# dt = min (dtAdv, dtCfl, dtStretch), where dtAdv is equal to LCFL / |gradU|
-# For dtAdv, the possible choices are the following:
-# 'vort' (infinite norm of vorticity) : dtAdv = LCFL / |vort|
-# 'gradU' (infinite norm of velocity gradient), dtAdv = LCFL / |gradU|
-# 'deform' (infinite norm of deformation tensor),
-# dtAdv = LCFL / (0.5(gradU + gradU^T))
-# The operator is defined on a Cartesian topology with ghost points, that will
-# be dedicated to operators using finite differences.
-from parmepy.operator.adapt_timestep import AdaptTimeStep
-from parmepy.mpi.topology import Cartesian
-NBGHOSTS = 2
-ghosts = np.ones((box.dimension)) * NBGHOSTS
-topostr = Cartesian(box, box.dimension, nbElem, ghosts=ghosts)
-
-from parmepy.operator.advection import Advection
-op['advection'] = Advection(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- method={Scales: 'p_M6',
- Splitting: 'classic'}) # Scales advection
-
-op['dtAdapt'] = AdaptTimeStep(velo, vorti,
- dt_adapt=dt,
- method={TimeIntegrator: RK3,
- SpaceDiscretisation: FD_C_4,
- dtCrit: ['vort', 'stretch']},
- topo=topostr,
- io_params={},
- lcfl=0.125,
- cfl=0.5)
-
-from parmepy.operator.advection import Advection
-op['advection'] = Advection(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- method={Scales: 'p_M6',
- Splitting: 'classic'}) # Scales advection
-
-from parmepy.operator.stretching import Stretching
-op['stretching'] = Stretching(velo, vorti, topo=topostr)
-
-from parmepy.operator.diffusion import Diffusion
-op['diffusion'] = Diffusion(vorti, resolutions={vorti:nbElem}, viscosity=VISCOSITY)
-
-from parmepy.operator.poisson import Poisson
-op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem})
-
-from parmepy.operator.differential import Curl
-op['curl'] = Curl(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
- method={SpaceDiscretisation: FD_C_4})
-
-# Discretisation of computational operators and link
-# to the associated topologies.
-for ope in op.values():
- ope.discretize()
-topofft = op['poisson'].discreteFields[op['poisson'].vorticity].topology
-topocurl = op['curl'].discreteFields[op['curl'].invar].topology
-topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
-topostr = op['stretching'].discreteFields[op['stretching'].velocity].topology
-
-# Penalization of the velocity on a sphere inside the domain.
-# We enforce penalization on the same topology as for fft operators.
-from parmepy.operator.penalization import Penalization
-op['penalization'] = Penalization(variables=velo, obstacles=[sphere], coeff=[1e8],
- topo=topofft)
-op['penalization'].discretize()
-
-orig = box.origin.copy()
-entr_length = box.length.copy()
-step = topofft.mesh.space_step
-entr_length[0] = 8 * step[0]
-entrance = ControlBox(domain=box, origin=orig, lengths=entr_length)
-op['killvorticity'] = Penalization(variables=vorti, obstacles=[entrance], coeff=[1e12],
- topo=topofft)
-op['killvorticity'].discretize()
-
-# Correction of the velocity, used as post-process for Poisson solver,
-# based on a fixed flowrate through the entrance of the domain.
-# Time-dependant required-flowrate (Variable Parameter)
-def computeFlowrate(simu):
- # === Time-dependant flow rate ===
- t = simu.tk
- Tstart = 3.0
- flowrate = np.zeros(3)
- flowrate[0] = uinf * box.length[1] * box.length[2]
- if t >= Tstart and t <= Tstart + 1.0:
- flowrate[1] = sin(pi * (t - Tstart)) * \
- box.length[1] * box.length[2]
- # === Constant flow rate ===
- # flowrate = np.zeros(3)
- # flowrate[0] = uinf * box.length[1] * box.length[2]
- return flowrate
-req_flowrate = VariableParameter(formula=computeFlowrate)
-
-frate = npw.zeros(3)
-frate[0] = uinf * box.length[1] * box.length[2]
-req_flowrate = VariableParameter(data=frate, name='flowRate')
-
-from parmepy.operator.velocity_correction import VelocityCorrection
-op['correction'] = VelocityCorrection(velo, vorti,
- req_flowrate=req_flowrate, topo=topofft)
-op['correction'].discretize()
-
-# ========= Bridges between the different topologies =========
-from parmepy.operator.redistribute import Redistribute
-distr = {}
-distr['fft2curl'] = Redistribute(variables=[velo], opFrom=op['penalization'], opTo=op['curl'])
-distr['fft2advec'] = Redistribute(variables=[velo, vorti], opFrom=op['poisson'], opTo=op['advection'])
-distr['curl2fft'] = Redistribute(op['curl'], op['poisson'], variables=[vorti])
-distr['adv2str'] = Redistribute(op['advection'], op['stretching'], variables=[velo, vorti])
-distr['str2diff'] = Redistribute(variables=[vorti], opFrom=op['stretching'], opTo=op['diffusion'])
-distr['curl2adv'] = Redistribute(variables=[velo, vorti], opFrom=op['curl'], opTo=op['advection'])
-distr['curl2str'] = Redistribute(variables=[velo, vorti], opFrom=op['curl'], opTo=op['stretching'])
-distr['fft2str'] = Redistribute(variables=[velo, vorti], opFrom=op['poisson'], opTo=op['stretching'])
-distr['str2curl'] = Redistribute(variables=[velo], opFrom=op['stretching'], opTo=op['curl'])
-for ope in distr.values():
- ope.discretize()
-
-# ========= Monitoring operators =========
-
-from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
- Reprojection_criterion
-monitors = {}
-monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5, prefix='fft_io')
-monitors['printerFFT0'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5, prefix='fft_io0')
-
-monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
- formattype=HDF5, prefix='curl_io')
-
-monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5, prefix='advec_io',
- xmfalways=True)
-
-monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
- io_params={},
- viscosity=VISCOSITY, isNormalized=False)
-
-reprojection_constant = 0.04
-reprojection_rate = 1
-monitors["reprojection"] = Reprojection_criterion(vorti, reprojection_constant,
- reprojection_rate, topo=topofft,
- io_params={})
-# Add reprojection for poisson operator
-op["poisson"].activateProjection(monitors["reprojection"])
-
-# Compute and save drag and lift on a
-# 3D Control box
-ref_step = topostr.mesh.space_step
-cbpos = boxorigin.copy()#
-cbpos[0] += 10 * ref_step[0]
-cblength = boxlength.copy()#
-cblength[0] -= 20 * ref_step[0]
-cb = ControlBox(domain=box, origin=cbpos, lengths=cblength)
-coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
-
-monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
- topo=topostr, volumeOfControl=cb,
- obstacles=[sphere], io_params={})
-
-step_dir = ref_step[0]
-thickSliceXY = ControlBox(domain=box, origin=[-2.0, -2.56, -2.0 * step_dir],
- lengths=[4., 5.12, 4.0 * step_dir])
-## #thickSliceXY = ControlBox(box, origin=[-2.0, -pi, -2.0 * step_dir],
-## # lengths=[8.0*pi, 2.0*pi, 4.0 * step_dir])
-monitors['printerSliceXY'] = Printer(variables=[velo, vorti], topo=topofft,
- frequency=1, formattype=HDF5, prefix='sliceXY',
- subset=thickSliceXY, xmfalways=True)
-
-## thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
-## lengths=[10.24, 4.0 * step_dir, 5.12])
-## monitors['printerSliceXZ'] = Printer(variables=[velo, vorti], topo=topofft,
-## frequency=100, formattype=HDF5, prefix='sliceXZ',
-## subset=thickSliceXZ, xmfalways=True)
-
-## subBox = ControlBox(box, origin=[-0.7, -1.25, -1.25], lengths=[7.0, 2.5, 2.5])
-## monitors['printerSubBox'] = Printer(variables=[velo, vorti], topo=topofft,
-## frequency=100, formattype=HDF5, prefix='subBox',
-## subset=subBox, xmfalways=True)
-
-# ========= Setup for all declared operators =========
-for ope in op.values():
- ope.setUp()
-for ope in distr.values():
- ope.setUp()
-# Set up for monitors
-for monit in monitors.values():
- monit.setUp()
-
-allops = dict(op.items() + distr.items() + monitors.items())
-
-# ========= Simulation setup =========
-#simu = Simulation(tinit=0.0, tend=15., timeStep=0.0125, iterMax=1000000)
-simu = Simulation(tinit=0.0, tend=75., timeStep=dt['dt'], iterMax=70)
-
-ind = sphere.discretize(topofft)
-vdfft = velo.discreteFields[topofft].data
-wdfft = vorti.discreteFields[topofft]
-
-## ========= Fields initialization =========
-# Initialisation, mode 1:
-# - compute velocity on topofft
-# - penalize
-# - compute vorticity with curl
-def initFields_mode1():
- # The velocity is initialized on the fftw topology
- # penalized and distributed on curl topo
- velo.initialize(topo=topofft)
- op['penalization'].apply(simu)
- distr['fft2curl'].apply(simu)
- distr['fft2curl'].wait()
- # Vorticity is initialized after penalization of velocity
- op['curl'].apply(simu)
- distr['curl2adv'].apply(simu)
- distr['curl2adv'].wait()
-
-## Call initialization
-initFields_mode1()
-
-##### Check initialize process results #####
-norm_vel_fft = velo.norm(topofft)
-norm_vel_curl = velo.norm(topocurl)
-norm_vort_fft = vorti.norm(topofft)
-norm_vort_curl = vorti.norm(topocurl)
-
-#assert np.allclose(norm_vort_curl, norm_vort_fft)
-assert np.allclose(norm_vel_curl, norm_vel_fft)
-
-# Print initial state
-for mon in monitors.values():
- mon.apply(simu)
-
-## Note Franck : tests init ok, same values on both topologies, for v and w.
-# === After init ===
-#
-# v init on topocurl (which may be equal to topofft)
-#
-# === Definition of the 'one-step' sequence of operators ===
-#
-# 1 - w = curl(v), topocurl
-# 2 - w = advection(v,w), topo-advec
-# 3 - w = stretching(v,w), topostr
-# 4 - w = diffusion(w), topofft
-# 5 - v = poisson(w), topofft
-# 6 - v = correction(v, w), topofft
-# 7 - v = penalization(v), topofft
-#
-# Required bridges :
-# 1 --> 2 : (v,w) on topo-advec
-# 2 --> 3 : (v,w) on topostr
-# 3 --> 4 : w on topofft
-# 7 --> 1 : v on topocurl
-
-#fullseq = ['advection', # in: v,w out: w, on topoadvec
-# 'adv2str', 'stretching', # in: v,w out: w on topostr
-# # in: w, out: v, w on topofft
-# 'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
-# 'penalization',
-# 'fft2curl', 'curl', 'forces',
-# 'curl2adv', 'energy', 'profiles', 'printerAdvec',
-# ]
-
-fullseq = ['stretching',
- 'str2diff', 'diffusion', 'poisson', 'correction',
- 'penalization',
- 'fft2curl', 'curl',
- 'curl2fft', 'poisson', 'correction',
- 'fft2adv',
- 'advection',
- #'printerSliceXY', #'printerSliceXZ', 'printerSubBox',
- #'energy',
- 'adv2str',
- 'forces',
- #'dtAdapt'
- ]
-
-cb2 = op['correction'].cb
-#cb.discretize(topo=topofft)
-def pseudo_norm(topo, var):
- sloc = npw.zeros(3)
- gloc = npw.zeros(3)
- sl = topo.mesh.iCompute
- for i in xrange(3):
- sloc[i] = np.sum((var[i][sl]))
- # sfft[i] = la.norm((vfft[i][slfft]))
- topo.comm.Allreduce(sloc, gloc)
- return gloc
-
-def run(sequence):
- ## for name in sequence:
- ## assert allops.keys().count(name) > 0 or distr.keys().count(name) > 0\
- ## or monitors.keys().count(name) > 0, 'unknow key:' + name
- ## allops[name].apply(simu)
- op['advection'].apply(simu)
- #monitors['printerFFT0'].apply(simu)
- op['killvorticity'].apply(simu)
- #monitors['printerFFT'].apply(simu)
-
- distr['adv2str'].apply(simu)
- distr['adv2str'].wait()
- op['stretching'].apply(simu)
- distr['str2diff'].apply(simu)
- distr['str2diff'].wait()
- op['diffusion'].apply(simu)
- op['poisson'].apply(simu)
- ## tata=np.zeros(3)
- ## titi=np.zeros(3)
- ## toto = np.zeros(3)
- ## #tutu = np.zeros(3)
- ## sl = cb2.slices[topofft]
- ## for dime in xrange(3):
- ## toto[dime] = np.sum(wdfft.data[dime][sl])# wdfft.integrate_on_proc(cb2, component=dime)
- ## #tutu[dime] = wdfft.integrate_on_proc(cb2, useSlice=False,
- ## # component=dime)
- ## wdfft.topology.comm.Allreduce(toto, tata)
- ## # wdfft.topology.comm.Allreduce(toto, titi)
-
- ## print tata, titi
- ## print pseudo_norm(topofft, wdfft)
- op['correction'].apply(simu)
- op['penalization'].apply(simu)
- distr['fft2curl'].apply(simu)
- distr['fft2curl'].wait()
- op['curl'].apply(simu)
- distr['curl2fft'].apply(simu)
- distr['curl2fft'].wait()
- op['poisson'].apply(simu)
- op['correction'].apply(simu)
-
- distr['fft2advec'].apply(simu)
- distr['fft2str'].apply(simu)
- distr['fft2str'].wait()
- monitors['forces'].apply(simu)
- distr['fft2advec'].wait()
-
-
-seq = fullseq
-
-simu.initialize()
-
-def checkCorrec():
- vcor = op["correction"].discreteOperator
- cbtest = vcor.cb
- sref = cbtest.lowerS[0]
- integ = npw.zeros(6)
- globalres = npw.zeros(6)
- for i in xrange(3):
- integ[i] = vcor.velocity.integrateOnSurf_proc(sref, component=i)
- for i in xrange(3):
- integ[i + 3] = vcor.vorticity.integrate_on_proc(cbtest, component=i)
- vcor.velocity.topology.comm.Allreduce(integ, globalres)
- print vcor.topo.rank, 'integ ...', integ
- print vcor.topo.rank, 'integ global ...', integ
-
-while not simu.isOver:
- if topocurl.rank == 0:
- simu.printState()
- run(seq)
- checkCorrec()
- simu.advance()
-
-
-
-## print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
-
-# === Finalize for all declared operators ===
-
-# Note FP : bug in fft finalize. To be checked ...
-# --> clean_fftw must be called only once
-#for ope in op.values():
-# ope.finalize()
-## Clean memory buffers
-
-fftw2py.clean_fftw_solver(box.dimension)
-for ope in distr.values():
- ope.finalize()
-for monit in monitors.values():
- monit.finalize()
diff --git a/Examples/TaylorGreen3D_debug.py b/Examples/TaylorGreen3D_debug.py
new file mode 100644
index 0000000000000000000000000000000000000000..80578cf57f1a0e81a9246bf195120f269e131ef8
--- /dev/null
+++ b/Examples/TaylorGreen3D_debug.py
@@ -0,0 +1,276 @@
+#!/usr/bin/python
+
+"""
+Taylor Green 3D : see paper van Rees 2011.
+
+All parameters are set and defined in python module dataTG.
+
+"""
+
+import parmepy as pp
+from parmepy.f2py import fftw2py
+import numpy as np
+from parmepy.fields.continuous import Field
+from parmepy.fields.variable_parameter import VariableParameter
+from parmepy.operator.advection import Advection
+from parmepy.operator.stretching import Stretching
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.operator.adapt_timestep import AdaptTimeStep
+from parmepy.operator.redistribute_intra import RedistributeIntra
+from parmepy.operator.hdf_io import HDF_Writer
+from parmepy.operator.energy_enstrophy import EnergyEnstrophy
+from parmepy.problem.simulation import Simulation
+from parmepy.methods_keys import Scales, TimeIntegrator,\
+ Splitting, dtCrit, SpaceDiscretisation
+from parmepy.numerics.integrators.runge_kutta3 import RK3 as RK3
+from parmepy.numerics.finite_differences import FD_C_4
+from parmepy.mpi import main_rank, MPI
+from parmepy.tools.parameters import Discretization, IO_params
+
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+
+# ====== pi constant and trigonometric functions ======
+pi = np.pi
+cos = np.cos
+sin = np.sin
+
+VISCOSITY = 1. / 1600.
+# ======= Domain =======
+dim = 3
+Nx = Ny = Nz = 65
+g = 2
+box = pp.Box(length=[2.0 * pi, 2.0 * pi, 2.0 * pi])
+
+# A global discretization with ghost points
+d3dg = Discretization([Nx, Ny, Nz], [g, g, g])
+# A global discretization, without ghost points
+d3d = Discretization([Nx, Ny, Nz])
+# Default topology (i.e. 3D, with ghosts)
+topo_with_ghosts = box.create_topology(d3dg)
+
+# ======= Function to compute TG velocity =======
+def computeVel(res, x, y, z, t):
+ res[0][...] = sin(x) * cos(y) * cos(z)
+ res[1][...] = - cos(x) * sin(y) * cos(z)
+ res[2][...] = 0.
+ return res
+
+# ======= Function to compute reference vorticity =======
+def computeVort(res, x, y, z, t):
+ res[0][...] = - cos(x) * sin(y) * sin(z)
+ res[1][...] = - sin(x) * cos(y) * sin(z)
+ res[2][...] = 2. * sin(x) * sin(y) * cos(z)
+ return res
+
+# ======= Fields =======
+velo = Field(domain=box, formula=computeVel,
+ name='Velocity', isVector=True)
+vorti = Field(domain=box, formula=computeVort,
+ name='Vorticity', isVector=True)
+
+# ======= Parameter Variable (adaptative timestep) =======
+dt = VariableParameter(data=0.0125, name='dt')
+
+# ========= Simulation setup =========
+simu = Simulation(tinit=0.0, tend=10.0, timeStep=0.0125, iterMax=50)
+
+# Adaptative timestep method : dt = min(values(dtCrit))
+# where dtCrit is a list of criterions on which the computation
+# of the adaptative time step is based
+# ex : dtCrit = ['gradU', 'cfl', 'stretch'], means :
+# dt = min (dtAdv, dtCfl, dtStretch), where dtAdv is equal to LCFL / |gradU|
+# For dtAdv, the possible choices are the following:
+# 'vort' (infinite norm of vorticity) : dtAdv = LCFL / |vort|
+# 'gradU' (infinite norm of velocity gradient), dtAdv = LCFL / |gradU|
+# 'deform' (infinite norm of deformation tensor),
+# dtAdv = LCFL / (0.5(gradU + gradU^T))
+op = {}
+iop = IO_params("time_step")
+
+op['dtAdapt'] = AdaptTimeStep(velo, vorti, simulation=simu,
+ discretization=topo_with_ghosts,
+ method={TimeIntegrator: RK3,
+ SpaceDiscretisation: FD_C_4,
+ dtCrit: ['deform', 'cfl', 'stretch']},
+ io_params=iop,
+ lcfl=0.125,
+ cfl=0.5)
+
+op['advection'] = Advection(velo, vorti,
+ discretization=d3d,
+ method={Scales: 'p_M6',
+ Splitting: 'classic'}
+ )
+
+op['stretching'] = Stretching(velo, vorti,
+ discretization=topo_with_ghosts)
+
+op['diffusion'] = Diffusion(viscosity=VISCOSITY, vorticity=vorti,
+ discretization=d3d)
+
+op['poisson'] = Poisson(velo, vorti, discretization=d3d, projection=1)
+
+# ===== Discretization of computational operators ======
+for ope in op.values():
+ ope.discretize()
+
+topofft = op['poisson'].discreteFields[vorti].topology
+topoadvec = op['advection'].discreteFields[vorti].topology
+
+# ==== Operators to map data between the different computational operators ===
+# (i.e. between topologies)
+distr = {}
+distr['adv2str'] = RedistributeIntra(source=op['advection'],
+ target=op['stretching'],
+ variables=[velo, vorti])
+distr['str2adv'] = RedistributeIntra(source=op['stretching'],
+ target=op['advection'],
+ variables=[velo, vorti])
+distr['fft2str'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretching'],
+ variables=[velo])
+distr['fft2str2'] = RedistributeIntra(source=op['poisson'],
+ target=op['stretching'],
+ variables=[velo, vorti])
+distr['str2fft'] = RedistributeIntra(source=op['stretching'],
+ target=op['diffusion'],
+ variables=[vorti])
+# ## ========= Monitoring operators =========
+
+iop = IO_params('TG_io', frequency=100)
+writer = HDF_Writer(variables={velo: topofft, vorti: topofft},
+ io_params=iop)
+
+io_ener=IO_params('energy_enstrophy')
+energy = EnergyEnstrophy(velo, vorti, discretization=topofft,
+ io_params=io_ener)
+writer.discretize()
+energy.discretize()
+
+# ========= Setup for all declared operators =========
+time_setup = MPI.Wtime()
+for ope in op.values():
+ ope.setup()
+for ope in distr.values():
+ ope.setup()
+writer.setup()
+energy.setup()
+
+print '[', main_rank, '] total time for setup:', MPI.Wtime() - time_setup
+
+# ========= Fields initialization =========
+# - initialize velo + vort on topostr
+# - penalize vorticity
+# - redistribute topostr --> topofft
+
+time_init = MPI.Wtime()
+
+
+def initFields_mode1():
+ velo.initialize(topo=topoadvec)
+ vorti.initialize(topo=topoadvec)
+
+
+initFields_mode1()
+print '[', main_rank, '] total time for init :', MPI.Wtime() - time_init
+
+fullseq = []
+
+
+def run(sequence):
+ op['advection'].apply(simu)
+ distr['adv2str'].apply(simu)
+ distr['adv2str'].wait()
+ op['stretching'].apply(simu)
+ distr['str2fft'].apply(simu)
+ distr['str2fft'].wait()
+ op['diffusion'].apply(simu)
+ op['poisson'].apply(simu)
+ energy.apply(simu)
+ writer.apply(simu)
+ distr['fft2str2'].apply(simu)
+ distr['fft2str2'].wait()
+ op['dtAdapt'].apply(simu)
+ op['dtAdapt'].wait()
+ distr['str2adv'].apply(simu)
+ distr['str2adv'].wait()
+
+# ==== Serialize some data of the problem to a "restart" file ====
+# def dump(filename):
+# """
+# Serialize some data of the problem to file
+# (only data required for a proper restart, namely fields in self.input
+# and simulation).
+# @param filename : prefix for output file. Real name = filename_rk_N,
+# N being current process number. If None use default value from problem
+# parameters (self.filename)
+# """
+# if filename is not None:
+# filename = filename
+# filedump = filename + '_rk_' + str(main_rank)
+# db = parmesPickle(filedump, mode='store')
+# db.dump(simu, 'simulation')
+# velo.dump(filename, mode='append')
+# vorti.dump(filename, mode='append')
+
+# ## ====== Load some data of the problem from a "restart" file ======
+# def restart(filename):
+# """
+# Load serialized data to restart from a previous state.
+# self.input variables and simulation are loaded.
+# @param filename : prefix for downloaded file.
+# Real name = filename_rk_N, N being current process number.
+# If None use default value from problem
+# parameters (self.filename)
+# """
+# if filename is not None:
+# filename = filename
+# filedump = filename + '_rk_' + str(main_rank)
+# db = parmesPickle(filedump, mode='load')
+# simu = db.load('simulation')[0]
+# simu.start = simu.time - simu.timeStep
+# ite = simu.currentIteration
+# simu.initialize()
+# simu.currentIteration = ite
+# print 'simu', simu
+# print ("load ...", filename)
+# velo.load(filename)
+# vorti.load(filename)
+# return simu
+
+seq = fullseq
+
+simu.initialize()
+#doDump = False
+#doRestart = False
+#dumpFreq = 5000
+#io_default={"filename":'restart'}
+#dump_filename = io.Writer(params=io_default).filename
+#===== Restart (if needed) =====
+# if doRestart:
+# simu = restart(dump_filename)
+# # Set up for monitors and redistribute
+# for ope in distr.values():
+# ope.setUp()
+# for monit in monitors.values():
+# monit.setUp()
+
+# ======= Time loop =======
+time_run = MPI.Wtime()
+while not simu.isOver:
+ if topofft.rank == 0:
+ simu.printState()
+ run(seq)
+ simu.advance()
+# # testdump = simu.currentIteration % dumpFreq is 0
+# # if doDump and testdump:
+# # dump(dump_filename)
+print '[', main_rank, '] total time for run :', MPI.Wtime() - time_run
+
+# ======= Finalize =======
+fftw2py.clean_fftw_solver(box.dimension)
+for ope in distr.values():
+ ope.finalize()
+writer.finalize()
+energy.finalize()