diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index 02605888087219a25527f7d4921b3e2d6bbcca8a..b3b25b93975e11597a78064961cfcc8074f5a14c 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -10,25 +10,7 @@ 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.fields.continuous import Field
-from parmepy.mpi.topology import Cartesian
-from parmepy.fields.variable_parameter import VariableParameter
-from parmepy.domain.obstacle.sphere import Sphere
-from parmepy.operator.advection import Advection
-from parmepy.operator.stretching import Stretching
-from parmepy.operator.poisson import Poisson
-from parmepy.operator.velocity_correction import VelocityCorrection
-#from parmepy.operator.vorticity_correction import VorticityCorrection
-from parmepy.operator.diffusion import Diffusion
-from parmepy.operator.penalization import Penalization
-from parmepy.operator.differential import Curl
-from parmepy.operator.adapt_timestep import AdaptTimeStep
-from parmepy.operator.redistribute import Redistribute
-from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
- Reprojection_criterion
from parmepy.problem.simulation import Simulation
-from parmepy.constants import VTK, HDF5
-from parmepy.domain.obstacle.planes import PlaneBoundaries, SubPlane
from parmepy.domain.obstacle.controlBox import ControlBox
from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
Remesh, Support, Splitting, dtCrit, SpaceDiscretisation, GhostUpdate
@@ -39,6 +21,7 @@ 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
@@ -50,44 +33,43 @@ pi = np.pi
cos = np.cos
sin = np.sin
-## Flow constants
+# ====== Flow constants =======
uinf = 1.0
VISCOSITY = 1. / 1600.
-## Domain
+# ========= Geometry of the domain =========
dim = 3
-#boxlength = npw.realarray([8.0 * pi, 2.0 * pi, 2.0 * pi])
-#boxorigin = npw.realarray([-2.0, -pi, -pi])
+#boxlength = npw.realarray([2.0, 2.0, 2.0])
+#boxorigin = npw.realarray([-1., -1., -1.])
boxlength = npw.realarray([10.24, 5.12, 5.12])
boxorigin = npw.realarray([-2.0, -2.56, -2.56])
+# The domain
box = pp.Box(dim, length=boxlength, origin=boxorigin)
-## Global resolution
+# Sphere inside the domain
+RADIUS = 0.5
+sphere_pos = [0., 0., 0.]
+from parmepy.domain.obstacle.sphere import Sphere
+sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
+
+# ========= Discretisation parameters =========
#nbElem = [1025, 513, 513]
nbElem = [129, 65, 65]
+#nb = 17
+#nbElem = [nb, nb, nb]
-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
+# ========= 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.
@@ -95,37 +77,16 @@ def computeVort(res, x, y, z, t):
res[2][...] = 0.
return res
-## Vector Fields
+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')
-## Parameter Variable (adaptative timestep)
-data = {'dt': 0.0125}
-dt_adapt = VariableParameter(data)
-
-## Usual Cartesian topology definition
-# At the moment we use two (or three?) topologies :
-# - "topo" for Stretching and all operators based on finite differences.
-# --> ghost layer = 2
-# - topo from Advection operator for all the other operators.
-# --> no ghost layer
-# - topo from fftw for Poisson and Diffusion.
-# Todo : check compat between scales and fft operators topologies.
-NBGHOSTS = 2
-ghosts = np.ones((box.dimension)) * NBGHOSTS
-topostr = Cartesian(box, box.dimension, nbElem,
- ghosts=ghosts)
-
-
-## === Obstacles (sphere + up and down plates) ===
-#sphere_pos = (boxlength - boxorigin) * 0.5
-RADIUS = 0.5
-sphere_pos = [0., 0., 0.]
-sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
+# ========= Operators that describe the problem =========
+op = {}
-## === Tools Operators ===
# 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
@@ -134,12 +95,20 @@ sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
# 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 = {}
+# '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)
+
op['dtAdapt'] = AdaptTimeStep(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- dt_adapt=dt_adapt,
+ dt_adapt=dt,
method={TimeIntegrator: RK3,
SpaceDiscretisation: FD_C_4,
dtCrit: ['vort', 'stretch']},
@@ -148,142 +117,151 @@ op['dtAdapt'] = AdaptTimeStep(velo, vorti,
lcfl=0.125,
cfl=0.5)
-## === Computational Operators ===
+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,
resolutions={velo: nbElem, vorti: nbElem},
topo=topostr)
+from parmepy.operator.diffusion import Diffusion
op['diffusion'] = Diffusion(vorti, resolution=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, GhostUpdate: True})
+ method={SpaceDiscretisation: FD_C_4})
-## Discretisation of computational operators
+# Discretisation of computational operators and link
+# to the associated topologies.
for ope in op.values():
ope.discretize()
-
-# Get fft topology and use it for penalization and correction operators
-opfft = op['poisson']
-topofft = opfft.discreteFields[opfft.vorticity].topology
+topofft = op['poisson'].discreteFields[op['poisson'].vorticity].topology
topocurl = op['curl'].discreteFields[op['curl'].invar].topology
topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
-# space step
-ref_step = topofft.mesh.space_step
-step_dir = ref_step[0]
+topostr = op['stretching'].discreteFields[op['stretching'].velocity].topology
-## Penalization velocity
+# 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(velo, [sphere], coeff=[1e8],
topo=topofft, resolutions={velo: nbElem})
+op['penalization'].discretize()
-## Time-dependant required-flowrate (Variable Parameter)
+# 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,
resolutions={velo: nbElem,
vorti: nbElem},
req_flowrate=req_flowrate, topo=topofft)
-
-op['penalization'].discretize()
op['correction'].discretize()
-## === Bridges between the different topologies ===
+# ========= Bridges between the different topologies =========
+from parmepy.operator.redistribute import Redistribute
distr = {}
-
distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
+distr['fft2advec'] = Redistribute([velo, vorti],
+ op['poisson'], op['advection'])
distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
distr['adv2str'] = Redistribute([velo, vorti],
- op['advection'], op['stretching'])
+ op['advection'], op['stretching'])
distr['str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
distr['curl2adv'] = Redistribute([velo, vorti], op['curl'], op['advection'])
distr['curl2str'] = Redistribute([velo, vorti], op['curl'], op['stretching'])
-
-#distr['adv2corr'] = Redistribute([velo, vorti],
-# op['advection'], op['correc_vorti'])
-
+distr['fft2str'] = Redistribute([velo, vorti], op['poisson'], op['stretching'])
+distr['str2curl'] = Redistribute([velo], op['stretching'], op['curl'])
for ope in distr.values():
ope.discretize()
-# === Simulation setup ===
-#simu = Simulation(tinit=0.0, tend=75., timeStep=0.0125, iterMax=1000000)
-simu = Simulation(tinit=0.0, tend=75., timeStep=dt_adapt['dt'], iterMax=1000000)
+# ========= Monitoring operators =========
-## === Monitoring operators ===
+from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
+ Reprojection_criterion
monitors = {}
-# Save velo and vorti values on topofft
monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5)
+ formattype=HDF5, prefix='fft_io')
monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
- formattype=HDF5, prefix='curl_io_vorti')
+ formattype=HDF5, prefix='curl_io')
monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
formattype=HDF5, prefix='advec_io',
xmfalways=True)
-thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir],
- lengths=[10.24, 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=20, 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)
-# Compute and save enstrophy/Energy, on topofft
-io_default = {}
monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
io_params={},
viscosity=VISCOSITY, isNormalized=False)
-# Compute velocity and vorticity magnitude on topofft and print it in data file
-#io_default = {}
-#monitors['profiles'] = ProfilesOutput(velo, vorti, obstacles=[sphere],
-# topo=topofft, io_params={"frequency":1000})
-
-# Setup for reprojection
reprojection_constant = 0.04
reprojection_rate = 1
monitors["reprojection"] = Reprojection_criterion(vorti, reprojection_constant,
reprojection_rate, topofft,
io_params={})
-
# Add reprojection for poisson operator
op["poisson"].activateProjection(monitors["reprojection"])
-ind = sphere.discretize(topofft)
-vdfft = velo.discreteFields[topofft].data
-
-# Compute and save drag and lift on topofft
+# Compute and save drag and lift on a
# 3D Control box
-cbpos = boxorigin + 10 * ref_step
-cblength = boxlength - 25 * ref_step
+ref_step = topostr.mesh.space_step
+cbpos = boxorigin#
+cbpos[0] += 10 * ref_step[0]
+cblength = boxlength#
+cblength[0] -= 20 * ref_step[0]
cb = ControlBox(box, cbpos, cblength)
coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
-#print coeffForce
-# Monitor for forces
monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
- topostr, cb, obstacles=[sphere], io_params={})
-#monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-# topo, cb, io_params={})
+ topostr, cb, obstacles=[sphere], io_params={})
-# === Setup for all declared operators ===
+step_dir = ref_step[0]
+## thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir],
+## lengths=[10.24, 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=20, 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():
@@ -294,8 +272,15 @@ for monit in monitors.values():
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=5)
-## === Fields initialization ===
+ind = sphere.discretize(topofft)
+vdfft = velo.discreteFields[topofft].data
+wdfft = vorti.discreteFields[topofft]
+
+## ========= Fields initialization =========
# Initialisation, mode 1:
# - compute velocity on topofft
# - penalize
@@ -306,14 +291,11 @@ def initFields_mode1():
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)
- op['advection'].apply(simu)
- # Distribute vorticity on topostr
- distr['adv2str'].apply(simu)
- # From this point both velocity and vorticity are initialized
- # on topocurl and topoadvec. velocity is also init. on topofft.
+ distr['curl2adv'].wait()
## Call initialization
initFields_mode1()
@@ -328,8 +310,8 @@ norm_vort_curl = vorti.norm(topocurl)
assert np.allclose(norm_vel_curl, norm_vel_fft)
# Print initial state
-#for mon in monitors.values():
-# mon.apply(simu)
+for mon in monitors.values():
+ mon.apply(simu)
## Note Franck : tests init ok, same values on both topologies, for v and w.
# === After init ===
@@ -366,17 +348,68 @@ fullseq = ['stretching',
'penalization',
'fft2curl', 'curl',
'curl2fft', 'poisson', 'correction',
+ 'fft2adv',
'advection',
- 'printerSliceXY', #'printerSliceXZ', 'printerSubBox',
- 'energy', 'adv2str', 'forces',
- 'dtAdapt']
-
-
+ #'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)
+ ## 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)
+ 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
@@ -398,6 +431,7 @@ while not simu.isOver:
#for ope in op.values():
# ope.finalize()
## Clean memory buffers
+
fftw2py.clean_fftw_solver(box.dimension)
for ope in distr.values():
ope.finalize()