Add new example

Examples/NS_planeJet.py

+#!/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.constants import ORDER, np, PARMES_REAL
+from parmepy.f2py import fftw2py
+import numpy as np
+from parmepy.fields.continuous import Field
+from parmepy.variables.variables import Variables
+from parmepy.mpi.topology import Cartesian
+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 import Redistribute
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.printer import Printer
+from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from parmepy.problem.simulation import Simulation
+from parmepy.operator.adapt_timestep import AdaptTimeStep
+from parmepy.methods_keys import Scales
+from parmepy.operator.differential import Curl
+from parmepy.numerics.updateGhosts import UpdateGhosts
+from parmepy.mpi.main_var import main_size
+from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
+    Remesh, Support, Splitting, dtAdvecCrit, 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 L4_2 as rmsh
+from parmepy.f2py import fftw2py
+
+# problem dimension
+dim = 3
+# resolution
+nb = 129
+# number of ghosts in usual cartesian topo
+NBGHOSTS = 2
+ADVECTION_METHOD = {Scales: 'p_66'}
+# ADVECTION_METHOD = {TimeIntegrator: RK2,
+#                     Interpolation: Linear,
+#                     Remesh: rmsh,
+#                     Support: '',
+#                     Splitting: 'o2_FullHalf'}
+VISCOSITY = 5e-6
+
+
+## ----------- A 3d problem -----------
+print " ========= Start Navier-Stokes 3D (Taylor Green benchmark) ========="
+## pi constant
+pi = np.pi
+cos = np.cos
+sin = np.sin
+exp = np.exp
+abs = np.abs
+tanh = np.tanh
+## Domain
+box = pp.Box(dim, length=[1., ]*3)
+
+## Global resolution
+nbElem = [nb] * dim
+randX = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+                   dtype=PARMES_REAL, order=ORDER) - 0.5
+randY = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+                   dtype=PARMES_REAL, order=ORDER) - 0.5
+randZ = np.asarray(np.random.random([nb - 1, nb - 1,(nb - 1)/main_size]),
+                   dtype=PARMES_REAL, order=ORDER) - 0.5
+
+# ##initjet.f
+# width = 0.01
+# ampl3 = 0.3
+# ampl2 = 0.
+# ampl = 0.05
+# def computeVel(res, x, y, z, t):
+#     yy = abs(y - 0.5)
+#     aux=(0.1-2.*yy)/(4.*width)
+#     res[0][...] =0.5*(1.+tanh(aux))*(1.+ampl3*sin(8.*pi*(x)))*(1.+ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randX)
+#     res[1][...] = ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randY
+#     res[2][...] = ampl*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randZ
+#     return res
+
+
+# def initScal(res, x, y, z, t):
+#     yy = abs(y - 0.5)
+#     aux=(0.1-2.*yy)/(4.*width)
+#     res[0][...] = 0.5*(1.+tanh(aux))*(1.+ampl2*exp(-abs((0.1-2.*yy)/(4.*width)**2))*randZ)
+#     return res
+
+
+# ## JCP initial condition
+def computeVel(res, x, y, z, t):
+    res[0][...] = 0.5*(1.-tanh((abs(y-0.5)-0.1/2.)/0.02))*(1.+0.3*sin(8*pi*x))
+    res[1][...] = 0.
+    res[2][...] = 0.
+    return res
+
+
+def initScal(res, x, y, z, t):
+    res[0][...] = 0.5*(1.-tanh((abs(y-0.5)-0.1/2.)/0.02))*(1.+0.3*sin(8*pi*x))
+    return res
+
+
+## Fields
+velo = Field(domain=box, formula=computeVel,
+             name='Velocity', isVector=True)
+vorti = Field(domain=box,
+              name='Vorticity', isVector=True)
+scal = Field(domain=box, formula=initScal,
+              name='PassiveScalar', isVector=False)
+
+## Variables
+dt_adapt = Variables(domain=box, name='adaptative time step',
+                     data=[0.01])
+
+## Variables
+## 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.
+ghosts = [NBGHOSTS] * box.dimension
+topo = Cartesian(box, box.dimension, nbElem,
+                 ghosts=ghosts)
+
+## Navier Stokes Operators
+advec = Advection(velo, vorti,
+                  resolutions={velo: nbElem,
+                               vorti: nbElem},
+                  method=ADVECTION_METHOD
+                  )
+advecScal = Advection(velo, scal,
+                      resolutions={velo: nbElem,
+                                   scal: nbElem},
+                      method=ADVECTION_METHOD
+                      )
+
+stretch = Stretching(velo, vorti,
+                     resolutions={velo: nbElem,
+                                  vorti: nbElem},
+                     topo=topo
+                     )
+
+diffusion = Diffusion(vorti,
+                      resolution=nbElem,
+                      viscosity=VISCOSITY
+                      )
+
+poisson = Poisson(velo, vorti,
+                  resolutions={velo: nbElem,
+                               vorti: nbElem},
+                  projection=[True, 1, False],
+                  )
+
+c = Curl(velo, vorti,
+         resolutions={velo: nbElem,
+                      vorti: nbElem},
+         method={SpaceDiscretisation: fftw2py, GhostUpdate: False},
+         )
+
+## Tools Operators
+dtAdapt = AdaptTimeStep(velo, vorti,
+                        resolutions={velo: nbElem,
+                                     vorti: nbElem},
+                        dt_adapt=dt_adapt,
+                        method={TimeIntegrator: RK3,
+                                SpaceDiscretisation: FD_C_4,
+                                dtAdvecCrit: 'deform'},
+                        topo=topo,
+                        lcfl=0.125,
+                        cfl=None, # 0.5,
+                        prefix='./res_PS/dt_2p')
+
+## Simulation
+simu = Simulation(tinit=0.0,
+                  tend=4.5, timeStep=dt_adapt,
+                  iterMax=10000)
+
+
+
+# Bridges between the different topologies in order to
+# redistribute data.
+# 1 -Advection to stretching
+toGhost_vorti = Redistribute([vorti], advec, stretch)
+toGhost_velo = Redistribute([velo], poisson, stretch)
+# 2 - Stretching to Poisson/Diffusion
+fromGhost_vorti = Redistribute([vorti], stretch, diffusion)
+# 3 - Poisson to TimeStep
+distrPoissTimeStep = Redistribute([velo, vorti], poisson, dtAdapt)
+
+
+#  Define the problem to solve
+pb = NSProblem(operators=[toGhost_velo, advecScal, advec,
+                          toGhost_vorti,
+                          stretch,
+                          fromGhost_vorti,
+                          diffusion,
+                          poisson,
+                          distrPoissTimeStep, dtAdapt],
+               simulation=simu, dumpFreq=-1)
+
+## Setting solver to Problem (only operators for computational tasks)
+pb.pre_setUp()
+## Diagnostics related to the problem
+topofft = poisson.discreteFields[poisson.vorticity].topology
+# energy = Energy_enstrophy(velo, vorti, isNormalized=True,
+#                           topo=topofft,
+#                           viscosity=VISCOSITY,
+#                           frequency=1,
+#                           prefix='./res_PS/energy_S128_2p.dat')
+
+# initialisation
+vorti.setTopoInit(topofft)
+velo.setTopoInit(topofft)
+scal.setTopoInit(topofft)
+#pb.addMonitors([energy])
+pb.setUp()
+
+c.discretize()
+c.setUp()
+c.apply(simu)
+
+# p = Printer(variables=[scal,],
+#             topo=topofft,
+#             frequency=1,
+#             prefix='./res_PS/scal_S128_2p',
+#             ext=".vtk")
+# pb.addMonitors([p])
+# p._printStep()
+
+# pf = Printer(variables=[velo, vorti],
+#             topo=topofft,
+#             frequency=1,
+#             prefix='./res_PS/fields_S128_2p',
+#             ext=".vtk")
+# pb.addMonitors([pf])
+# pf._printStep()
+
+
+def run():
+    pb.solve()
+
+## Solve problem
+from parmepy.mpi import MPI
+print "Start computation ..."
+time = MPI.Wtime()
+run()
+print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
+
+
+pb.finalize()
+## Clean memory buffers
+#fftw2py.clean_fftw_solver(box.dimension)