Skip to content
Snippets Groups Projects
Commit 1115d596 authored by Chloe Mimeau's avatar Chloe Mimeau
Browse files

update NSDebug benchmark

parent 00a743dc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Examples/NSDebug.py
+ 63
− 23
View file @ 1115d596
......@@ -22,6 +22,7 @@ from parmepy.operator.velocity_correction import VelocityCorrection
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
......@@ -29,11 +30,18 @@ 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 dataNS_bb import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
OUTPUT_FREQ, FILENAME, PROJ, LCFL, CFL, CURL_METHOD, \
TIMESTEP_METHOD, OBST_Ox, OBST_Oy, OBST_Oz, RADIUS
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
import parmepy.tools.io_utils as io
## ----------- A 3d problem -----------
print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
......@@ -42,25 +50,24 @@ pi = np.pi
cos = np.cos
sin = np.sin
## Flow constants
uinf = 1.0
VISCOSITY = 1. / 1600.
## 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([10.24, 5.12, 5.12])
boxorigin = npw.realarray([-2.0, -2.56, -2.56])
box = pp.Box(dim, length=boxlength, origin=boxorigin)
nbElem = [129, 65, 65]
## Global resolutionsimu
## Global resolution
#nbElem = [1025, 513, 513]
# Initial upstream flow velocity
uinf = 1.0
nbElem = [129, 65, 65]
def computeFlowrate(simu):
##Time-dependant flow rate
# === Time-dependant flow rate ===
t = simu.tk
Tstart = 3.0
flowrate = np.zeros(3)
......@@ -68,7 +75,7 @@ def computeFlowrate(simu):
if t >= Tstart and t <= Tstart + 1.0:
flowrate[1] = sin(pi * (t - Tstart)) * \
box.length[1] * box.length[2]
##constant flow rate
# === Constant flow rate ===
# flowrate = np.zeros(3)
# flowrate[0] = uinf * box.length[1] * box.length[2]
return flowrate
......@@ -94,6 +101,10 @@ velo = Field(domain=box, formula=computeVel,
vorti = Field(domain=box, formula=computeVort,
name='Vorticity', isVector=True)
## 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.
......@@ -102,20 +113,47 @@ vorti = Field(domain=box, formula=computeVort,
# --> 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)
## === Computational Operators ===
## === 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
# 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 = {}
op['dtAdapt'] = AdaptTimeStep(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
dt_adapt=dt_adapt,
method={TimeIntegrator: RK3,
SpaceDiscretisation: FD_C_4,
dtCrit: ['vort', 'stretch']},
topo=topostr,
io_params={},
lcfl=0.125,
cfl=0.5)
## === Computational Operators ===
op['advection'] = Advection(velo, vorti,
resolutions={velo: nbElem, vorti: nbElem},
method=ADVECTION_METHOD)
resolutions={velo: nbElem,
vorti: nbElem},
method={Scales: 'p_M6',
Splitting: 'classic'}) # Scales advection
op['stretching'] = Stretching(velo, vorti,
resolutions={velo: nbElem, vorti: nbElem},
......@@ -126,7 +164,7 @@ op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem})
op['curl'] = Curl(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
method=CURL_METHOD)
method={SpaceDiscretisation: FD_C_4, GhostUpdate: True})
## Discretisation of computational operators
for ope in op.values():
......@@ -173,7 +211,8 @@ 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=0.0125, iterMax=1000000)
simu = Simulation(tinit=0.0, tend=75., timeStep=dt_adapt['dt'], iterMax=1000000)
## === Monitoring operators ===
monitors = {}
......@@ -193,7 +232,7 @@ thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.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=100, formattype=HDF5, prefix='sliceXY',
frequency=20, formattype=HDF5, prefix='sliceXY',
subset=thickSliceXY, xmfalways=True)
thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
......@@ -289,8 +328,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 ===
......@@ -328,8 +367,9 @@ fullseq = ['stretching',
'fft2curl', 'curl',
'curl2fft', 'poisson', 'correction',
'advection',
# 'printerSliceXY', 'printerSliceXZ', 'printerSubBox',
'energy', 'adv2str', 'forces']
'printerSliceXY', #'printerSliceXZ', 'printerSubBox',
'energy', 'adv2str', 'forces',
'dtAdapt']
def run(sequence):
......
Examples/dataNS_bb.py deleted 100644 → 0
+ 0
− 60
View file @ 00a743dc
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
from parmepy.f2py import fftw2py
import numpy as np
# problem dimension
dim = 3
# resolution
nb = 33
# number of ghosts in usual cartesian topo
NBGHOSTS = 2
# Adaptative timestep method
# For the "dtCrit" key, the possible choices are the following:
# 'vort' : def of dt_advection based on infinite norm of vorticity,
# 'gradU' : def of dt_advection based on infinite norm of velocity gradient,
# 'deform' : def of dt_advection based on infinite norm of deformation tensor
TIMESTEP_METHOD = {TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4,
dtCrit: 'deform'}
# Lagrangian CFL for adaptative timestep
LCFL = 0.125
# CFL (if CFL is None, no CFL condition is taken into account
# in the computation of adaptative timesteps)
CFL = 0.5
# Advection method
ADVECTION_METHOD = {Scales: 'p_M6', Splitting: 'classic'}
#ADVECTION_METHOD = {TimeIntegrator: RK2,
# Interpolation: Linear,
# Remesh: rmsh,
# Support: '',
# Splitting: 'o2_FullHalf'}
# Curl method
CURL_METHOD = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
#CURL_METHOD = {SpaceDiscretisation: 'fftw'}
# Viscosity of the fluid
VISCOSITY = 1. / 300.
# Vorticity projection (list of 2 elements):
# - the 1st element determines if the reprojection is needed
# - the 2nd element determines the reprojection frequency
# is terms of number of iterations IF the 1st element is True
# ex : PROJ = [True, 10] means that the reprojection
# of vorticty will be applied every 10 iterations
PROJ = [False, 10]
## pi constant
pi = np.pi
# Obstacle description
OBST_Ox = 0.
OBST_Oy = 0.
OBST_Oz = 0.
RADIUS = 0.5
# post treatment ...
OUTPUT_FREQ = 1
FILENAME = './res/energy_NS_sphere_FFTcurl_8cpu.dat'
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment