diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index 2047a42e188c7a1439f407db7dcd809dcc86e8bf..aeb96c8ca3d965af4c5cae2eab9d99e97f10b80c 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -41,12 +41,12 @@ cos = np.cos
sin = np.sin
## Domain
-boxlength = npw.realarray([2.0 * pi, pi, pi])
+boxlength = npw.realarray([2.0 * pi, 2 * pi, 2 * pi])
boxorigin = npw.realarray([0., 0., 0.])
box = pp.Box(dim, length=boxlength, origin=boxorigin)
-## Global resolution
-nbElem = [65, 33, 33]
+## Global resolutionsimu
+nbElem = [33,]*3
# Upstream flow velocity
uinf = 1.0
@@ -54,7 +54,7 @@ uinf = 1.0
# Function to compute potential flow past a sphere
def computeVel(res, x, y, z, t):
res[0][...] = uinf
- res[1][...] = 0.0
+ res[1][...] = 0.
res[2][...] = 0.
return res
@@ -115,6 +115,7 @@ for ope in op.values():
opfft = op['poisson']
topofft = opfft.discreteFields[opfft.vorticity].topology
topocurl = op['curl'].discreteFields[op['curl'].invar].topology
+topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
op['penalization'] = Penalization(velo, [sphere, bc], coeff=[1e8],
topo=topofft, resolutions={velo: nbElem})
@@ -147,12 +148,13 @@ distr['adv2str'] = Redistribute([velo, vorti],
# 3 - Stretching to Diffusion
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'])
for ope in distr.values():
ope.discretize()
# === Simulation setup ===
-simu = Simulation(tinit=0.0, tend=1, nbIter=10, iterMax=1000000)
+simu = Simulation(tinit=0.0, tend=10., nbIter=2000, iterMax=1000000)
## === Monitoring operators ===
monitors = {}
@@ -162,9 +164,13 @@ monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
formattype=HDF5, prefix='curl_io_vorti')
+monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
+ formattype=HDF5, prefix='advec_io',
+ xmfalways=True)
# Compute and save enstrophy/Energy, on topofft
io_default = {}
-monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft, io_params={},
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
+ io_params={},
viscosity=VISCOSITY, isNormalized=False)
# Setup for reprojection
reprojection_constant = 0.04
@@ -182,14 +188,18 @@ vdfft = velo.discreteFields[topofft].data
# Compute and save drag and lift on topofft
# 3D Control box
ref_step = topofft.mesh.space_step
-cbpos = boxorigin + boxlength + 10 * ref_step
+cbpos = boxorigin + 10 * ref_step
cblength = boxlength - 20 * ref_step
cb = ControlBox(box, cbpos, cblength)
-nu = 0.1
-coeffForce = 1.
+coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
+
+print coeffForce
# Monitor for forces
-monitors['forces'] = DragAndLift(velo, vorti, nu, coeffForce, topofft, cb,
- obstacles=[sphere], io_params={})
+## monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+## topo, cb,
+## obstacles=[sphere], io_params={})
+monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+ topo, cb, io_params={})
# === Setup for all declared operators ===
for ope in op.values():
@@ -217,11 +227,9 @@ def initFields_mode1():
# Vorticity is initialized after penalization of velocity
op['curl'].apply(simu)
# Distribute vorticity on topofft
- distr['curl2fft'].apply(simu)
+ distr['curl2str'].apply(simu)
# From this point both velocity and vorticity are initialized
- # on topofft and on topocurl
- # Remark : we only need to initialize v on topocurl
- # w computation and distribution is there just on test purpose.
+ # on topocurl and topoadvec. velocity is also init. on topofft.
## Call initialization
initFields_mode1()
@@ -237,11 +245,10 @@ print norm_vel_curl, norm_vel_fft
assert np.allclose(norm_vort_curl, norm_vort_fft)
assert np.allclose(norm_vel_curl, norm_vel_fft)
-## Print initial state
+# 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 ===
#
@@ -263,14 +270,23 @@ for mon in monitors.values():
# 3 --> 4 : w on topofft
# 7 --> 1 : v on topocurl
-fullseq = ['curl', 'printerCurl', # in: v, out : w both on topocurl
- 'curl2adv', 'advection', # in: v,w out: w, on topoadvec
- 'adv2str', 'stretching', # in: v,w out: w on topostr
+## 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',
+## 'curl2adv', 'energy', 'forces', 'printerAdvec',
+## ]
+fullseq = ['stretching', # in: v,w out: w on topostr
# in: w, out: v, w on topofft
'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
- 'penalization', 'energy', 'forces', 'printerFFT',
- 'fft2curl' # Back to topocurl for v
- ]
+ 'penalization',
+ 'fft2curl', 'curl',
+ 'curl2adv',
+ 'advection', # in: v,w out: w, on topoadvec
+ 'energy', 'printerAdvec',
+ 'adv2str', 'forces']
def run(sequence):
@@ -280,7 +296,9 @@ def run(sequence):
allops[name].apply(simu)
seq = fullseq
+
simu.initialize()
+
while not simu.isOver:
if topocurl.rank == 0:
simu.printState()
@@ -319,4 +337,3 @@ for ope in distr.values():
ope.finalize()
for monit in monitors.values():
monit.finalize()
-
diff --git a/Examples/dataNS_bb.py b/Examples/dataNS_bb.py
index e7783e24c32b401da71a320848754d8ff0f5af30..3652e6358dd1668bb82ab5d9150bb8f21e948669 100644
--- a/Examples/dataNS_bb.py
+++ b/Examples/dataNS_bb.py
@@ -36,7 +36,7 @@ ADVECTION_METHOD = {Scales: 'p_M6', Splitting: 'classic'}
# Splitting: 'o2_FullHalf'}
# Curl method
#CURL_METHOD = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
-CURL_METHOD = {SpaceDiscretisation: fftw2py, GhostUpdate: False}
+CURL_METHOD = {SpaceDiscretisation: fftw2py}
# Viscosity of the fluid
VISCOSITY = 1. / 300.
# Vorticity projection (list of 2 elements):