add a new example for fake 2D flow past a cylinder

Examples/NSDebug_faux2D.py

+#!/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. / 300.
+
+# ========= Geometry of the domain =========
+dim = 3
+boxlength = npw.realarray([6.0, 6.0, 0.08])
+boxorigin = npw.realarray([-2.0, -3.0, -0.04])
+#boxlength = npw.realarray([5.12, 5.12, 5.12])
+#boxorigin = npw.realarray([-2.56, -2.56, -2.56])
+# The domain
+box = pp.Box(dim, length=boxlength, origin=boxorigin)
+
+# Sphere inside the domain
+RADIUS = 0.5
+obst_pos = [0., 0., 0.]
+from parmepy.domain.obstacle.sphere import Sphere, HemiSphere
+from parmepy.domain.obstacle.cylinder import Cylinder, HemiCylinder
+#sphere = Sphere(box, position=obst_pos, radius=RADIUS)
+sphere = Cylinder(box, position=obst_pos, radius=RADIUS)
+
+# ========= Discretisation parameters =========
+nbElem = [301, 301, 5]
+#nbElem = [65, 65, 65]
+
+# ========= Vector Fields and variable parameters =========
+
+# Function to compute initial velocity
+def computeVel(res, x, y, z, t):
+    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)
+
+op['dtAdapt'] = AdaptTimeStep(velo, vorti,
+                              resolutions={velo: nbElem,
+                                           vorti: nbElem},
+                              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_M4', 
+                                    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})
+
+# 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(velo, [sphere], coeff=[1e8],
+                                  topo=topofft, resolutions={velo: nbElem})
+op['penalization'].discretize()
+
+# Kill the vorticity at the inlet.
+# We enforce penalization on the same topology as for fft operators.
+ref_step_fft = topofft.mesh.space_step
+ibpos = npw.zeros(dim)
+iblength = npw.zeros(dim)
+ibpos[...] = boxorigin[...]
+iblength[...] = boxlength[...]
+iblength[0] = 10 * ref_step_fft[0]
+inlet_band = ControlBox(box, ibpos, iblength)
+op['kill_vort'] = Penalization(vorti, [inlet_band], coeff=[1e10],
+                               topo=topofft, resolutions={vorti: nbElem})
+op['kill_vort'].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,
+                                      resolutions={velo: nbElem,
+                                                   vorti: nbElem},
+                                      req_flowrate=req_flowrate, topo=topofft)
+op['correction'].discretize()
+
+# ========= 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'])
+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['fft2str'] = Redistribute([velo, vorti], op['poisson'], op['stretching'])
+distr['str2curl'] = Redistribute([velo], op['stretching'], 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,
+                                 frequency=100, formattype=HDF5, prefix='fields',
+                                 xmfalways=True)
+
+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, 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 = npw.zeros(dim)
+cblength = npw.zeros(dim)
+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)
+
+monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+                                 topostr, cb, obstacles=[sphere], io_params={})
+
+step_dir = ref_step[0]
+#thickSliceXY = ControlBox(box, origin=[-2.56, -2.56, -2.0 * step_dir], 
+#                           lengths=[5.12, 5.12, 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=100., timeStep=0.0125, iterMax=1000000)
+#simu = Simulation(tinit=0.0, tend=75., timeStep=dt['dt'], iterMax=6)
+
+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)
+#    distr['adv2str'].apply(simu)
+#    distr['adv2str'].wait()    
+#    op['stretching'].apply(simu)
+#    distr['str2diff'].apply(simu)
+#    distr['str2diff'].wait()
+    op['diffusion'].apply(simu)
+    op['kill_vort'].apply(simu)
+    op['poisson'].apply(simu)
+    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)
+    monitors['printerFFT'].apply(simu)
+    monitors['energy'].apply(simu)
+    distr['fft2advec'].apply(simu)
+    distr['fft2str'].apply(simu)
+    distr['fft2str'].wait()
+    monitors['forces'].apply(simu)
+    distr['fft2advec'].wait()
+
+seq = fullseq
+
+simu.initialize()
+
+while not simu.isOver:
+    if topocurl.rank == 0:
+        simu.printState()
+    run(seq)
+    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()