Examples
my_first_polar.py
This example script my_first_polar.py is located in examples/. This code executes all different parts a polar computation could constitute. However you can of course break it up. Only launch computations in one script and then another that does the postprocessing, it is completely up to you.
import numpy as np
from e2dpolar.main import e2dpolar
from e2dpolar.in_out import gather_rundir_pkls, interp_bl_prof, load_dict, save_dict
import matplotlib.pyplot as plt
import sys
import os
# ======== STARTUP ============================================================
# Load instance
aero = e2dpolar()
# The object comes with a (huge) number of methods that allow you to
# mesh, run and postprocess CFD simulations. All inputs to the
# methods are specified through the object attributes. A full list
# can be found in the documentation or have a look at e2dpolar/main.py.
# Load the airfoil coordinates
aero.pf_coords = np.loadtxt('airfoils/naca63-418.dat')
# You do not need to use numpy here any 2D-array would do with x in the first
# and y in the second column. Nothing has happened yet, the coordinates were
# only loaded.
# ======== GRID ===============================================================
# Make surface mesh using PGL and volume grid with HypGrid2D.
# By default there are 512 points along the chord and 256 in the
# normal direction. The grid is partioned into blocks, by default
# the block size is 128. The total number of blocks in which the
# domain is partioned is 8. All this can be changed if wanted. The mesh
# really only needs to be generated once, so rerun=False would
# save some time and avoid you overwriting your old mesh.
# ----- Surface mesh ----
# The example airfoil coordinates used here do not have a closed
# trailing edge (TE), this is done by PGL automatically. Change this
# by setting aero.pf_te_thick_lim, which by default is set to
# None, ie nothing is done to the TE thickness. If the TE is not closed
# PGL will detect that and close it automatically. If you do set somehting
# here, then the original aerofoil coordindates will be altered
# (linear addition/removal of thickness from LE to TE) such that the TE
# thickness is 2e-3, which is standard for 1 m chord
# aerofoil in a wind tunnel. Let's say we want a specific TE thickness:
aero.pf_te_thick_lim = [2.00e-3, 'equal']
# By default the coordinates you provide will also be normalised by the
# chord (defined as the distance from LE to TE, where the LE is the point
# with max. distance from the TE) and rotated such that LE and TE are on
# the x-axis. This important as it changes the definition of the angle-
# of-attack. You can switch this off by setting the flags to False ie
# aero.pf_norm = False
# aero.pf_rotate_align = False
# by default these are set to True.
# Also note that by definition the coordinates should start at the TE
# and then move clockwise over the LE to the upper TE edge. If the
# direction is inverted, this should be caught and rectified by the
# program, but if not, the volume mesh generation will crash.
# If you want to use your own surface mesh set aero.af.points which
# holds the surface mesh from PGL. Just replace the first column by
# your x and the second by your y coordinates or use any of the PGL
# methods for altering surfaces meshes (also have a look at
# advanced_polars.py or step.py in examples/).
# You can control the inputs to the PGL AirfoilShape() class through
# aero.pf_pgl_input = {'even': False,
# 'dist': None,
# 'dLE': False,
# 'dTE': -1,
# 'close_te': True,
# 'linear': False}
# Checkout build/PGL/PGL/components/airfoil.py for all the options; 'dTE'
# might be interesting as it sets the number of grid points to be used
# to close the TE, if it is set to -1 as done by default, the TE grid
# spacing is set automatically similar to the spacing over the LE. Set
# this to an integer to set the exact no. of grid points to be used (uneven!).
# Let's make the surface mesh using PGL
rerun = True
aero.make_surface_mesh(rerun=rerun)
# The surface mesh and the AirfoilShape() used here from PGL are stored
# inside sfgrid/ together with a plot showing some of the mesh details and
# the surface gradients in 'surface_mesh.pdf'
# ----- Volume mesh -----
# HypGrid2D is used in growing the mesh from the surface and you can
# change the default inputs (they generally work quite well) and
# Basis2D is used to split the mesh into blocks. The volume grid and
# all associated files are located in the grid/ folder. You can open
# the grid.X2D file in ParaView or excute the 'pltgrid' command
# in the terminal to visualize the mesh using hpglview.
# aero.grid_bsize = 64 # would change the block size for example
aero.make_volume_mesh(rerun=rerun)
# ----- Surface + Volume
# There is a single command that generates the surface mesh using PGL
# and then grows a volume mesh from the surface using HypGrid2D:
# aero.make_mesh(rerun=True)
# ======== CFD SETUP ==========================================================
# Set the Reynolds number
aero.Re = 3e6
# All other flow parameters are derived from there, where the velocity,
# chord and density are assumed to be unity. You can change that but
# really shouldn't.
# Set the boundary layer flow case. For each a number of adequate
# default CFD run parameters is defined, all run with the k-omega SST model.
# Of course you can also change these default values. These are the three
# cases you may select from (only one can be active):
# aero.case = 'tr' # Transitional boundary layer, default with e^N + bypass
aero.case = 'tu' # Fully turbulent boundary layer
# through a roughness box.
# For larger Reynolds numbers (as early as 5e6) need to use the
# e^N + a bypass transition model (only if aero.case = 'tr'). You also got
# the option to use Langtry & Menter's gamma-Re_theta model.
# aero.tr_model = 'gamma-Retheta'.
# You should set the inflow turbulence level (only really relevant for
# transitional computations):
aero.ti = 0.001
# Finally set the angles-of-attack you would like to be computed
aero.aoas = np.arange(-4., 4., 2.)
# ======== Cluster
# Set the max computational time per aoa (hh:mm:ss). The computation is
# killed if it exceeds this.
aero.cpu_time = '00:10:00'
# Set the run queue, check availability beforehand
# aero.queue = 'workq'
aero.queue = 'windq'
# no. of processors (default: 1), if you like it faster increase the no. of
# processors but ensure that the no. of mesh blocks is divisble by the no. of
# procs. The blocks are given as:
# n_blocks = grid_ni/grid_bsize * grid_nj/grid_bsize
# , but are also outputted during the volume mesh generation.
aero.procs = 8
# ======== SIMULATIONS ========================================================
# In the directory where you have executed this script a directory will
# be created for the specified flow case Re<aero.Re>_<aero.case>, which then
# contains one CFD run directory per angle-of-attack a<aoa>. These run
# directories contain everything you need if you should want to rerun the
# computations and also keeps all the standard EllipSys2D outputs.
# Before launching the computations make sure to not drown the cluster in
# computations; do not make massive parameter sweeps and be aware of other
# people using the cluster (see documentation how to work on the cluster)!
# You can also change the convergence limit (reslim) for the different grid
# levels and number of iterations (mstep) to be used to converge the problem.
# The flow is solved at different grid levels to acclerate the convergence of
# the problem as you feed a good initial guess to the next more refined level.
# The finest grid is corresponds to level 1 and is equal to the one you specified
# each higher level reduces the no. of grid points by half. The default are 4
# grid levels, you will get output at each grid level, so you can easily check
# for grid sensitivity.
# The default should be fine, here it is changed to accelerate
# the computations on the testmachine
if aero._runmachine == 'testmachine':
aero.mstep = [10, 20, 30, 40]
aero.aoas = np.arange(-2., 6., 4.)
group_sims = False
else:
group_sims = True
aero.run_polar(group_sims=group_sims, rerun=rerun)
aero.check_finish()
# os.system('ls Re3.00e6_tu/a2.00')
# Wait for computations to finish. If you comment this then the code
# will continue and try to extract polars and do other postprocessing
# but will simply exit as the computations are not finished. If you set
# the rerun flags for meshing and for the CFD simulations to False, then
# the code will simply jump those tasks and continue with the post-
# processing. Try it!
# If the program complains later during data extraction or says the computations
# have not finished you should potentially increase the 'aero.cpu_time' to make
# sure that the computations fully converge/reach the desired residual level
# for all equations. Check the 'grid.res' file in the AoA folders to see the
# convergence history, maybe you can spot the issue. To help you the residuals
# and the convergence of the forces are plotted in pdfs, ie 'res.pdf', 'force_res.pdf'.
# You might have to tweak the simulation inputs to get it to converge
# or use boostconv by acdi!
# ======== POSTPROCESSING =====================================================
# ----- Extract polars -----
# you can extract the polar data you came for
aero.ext_polar(rerun=rerun)
# This method will also write a few plots to give an overview of the
# results. Additionally the polars are wirtten to grid-<grid_level>.clcd.
# Data is usually only extracted at the finest grid level (1) but you
# can also extract at different levels by setting aero.glevels = np.array([1,2]).
# Also note that there is much more data available from the computations,
# which are all stored in the object in the dictionary aero.polar. Surface data
# (like Cf and Cp) is stored for each angle-of-attack in aero.polar['surface_data']
# and the data ending up in the grid-<grid_level>.clcd in aero.polar['aero_data'].
# So you can use the data and output it as you like. The python directory
# holding the polar data is also stored in the run directory as 'polar.pkl'. So
# you can load it and do any kind of post-processing with it you like. Furthermore
# the polar and surface data is plotted and saved in figures inside the run
# directory in 'polars_gl<grid-level>.pdf' and 'profiles_gl<grid-level>.pdf', if
# running with transition you get also a pdf with the transition location.
# ----- Tools -----
# There are a multitude of methods available for post-processing, loading and
# gathering datasets as part of the pye2dpolar class, but also some that work
# independantly from it that are placed inside 'in_out.py' so check them out
# before making your own. One that might come handy when you have multiple
# run directories and you would like to assemble data for all Reynolds
# numbers or flow cases is a .pkl gathering method. By default it assembles
# the polar.pkl output files, but you can change that.
# Gather polars from run directories (Re*e6*), of course we only got one
# right now, but it is just to show how it works.
gather_rundir_pkls()
# ----- Extract flowfield -----
# If you like you can also extract the flow variables from the
# computations. Specify from which anlge-of-attack you would
# like to get the data. Be aware of memory contraints if you
# wish to extract a lot of fields and angles-of-attack. If you do not
# specify the fields all fields are extracted which should be there
# as a default. They are named as follows:
# ['u', 'v', 'p', 'vis', 'normal', 'vorticity', 'ko2_tke', 'ko2_omega']
# With the transition model on it is also 'tr1_gamma'.
# Here we only extrac the major ones the velocities and pressure,
# additionally we ask them to be plotted.
ext_fields = ['u', 'v', 'p']
# aero.ext_flowfield(-2.0, ext_fields=ext_fields, plt_fields=ext_fields,
# write_plot3d=False)
aero.ext_flowfield(2.0, ext_fields=ext_fields, plt_fields=ext_fields,
write_plot3d=False)
# You can also choose to get the flowfield as plot3d files, that can be
# opened with your favourite CFD data viewer. The data is
# again also saved in the object and can be found in the dict aero.flowfield
# including the vertex coordiantes of the grid x and y. You can loop over the
# routine as well. The dict is also saved in the AoA folder for which you
# extracted the flowfield inside 'flowfield.pkl' and contour plots are generated
# for the different fields and stored ('ff_<variable>_gl<grid-level>_aoa<AoA>.pdf').
# ----- Extract boundary layer -----
# There are in fact different ways you can extract boundary layers. One relies
# on having the flowfield data in memory or loaded as single block numpy array
# that you can provide to the method. So you need to have run aero.ext_flowfield()
# beforehand or need to load its output you have previously stored and provide it.
# Method #1
# Let's extract boundary layer profiles at four locations along the surface. The
# leading edge is at 0 and the x location is taken chord parallel with negative
# values on the pressure side.
xq = np.array([-0.4, -0.1, 0.2, 0.6])
# As we have the flowfield from AoA=2. in memory from running aero.ext_flowfield
# above we can simply run the method and get some profiles. This also generates
# some plots and output file inside the AoA folder to which the BL profiles
# belong. BLs are extracted normal to the surface and until where from the
# surface to extract the profiles is set by the user.
bl_profiles = aero.ext_bl(xq, output=True)
# This method is good if you want to quickly check BL profiles and compare
# between different AoAs or Reynolds numbers. Plots are stored under
# 'bl_profiles.pdf' in the AoA folder for which BLs where extracted.
# Method #2
# Another method that was developed by Franck Bertagnolio (frba) uses a
# modified version of the E2D postprocessor that extracts profiles all
# along the surface mesh and also computes BL characteristics along the
# aerofoil. Profiles are only extracted unitl the BL edge. As it
# extracts the BL profiles everywhere and needs to find the BL edge this
# method is slow (serial). BLs are extracted for all aoas so make sure
# you only extract those you need. It only works with the k-omega
# turbulence model right now.
# Let's just check the BLs of one AoA and also use some of the post-
# processing tools
aero.aoas = np.array([2.])
# If you have run this script before you might have the BL data already
# available so we can do some plotting, otherwise it will run the extraction
if not hasattr(aero, 'bl'):
try:
# this method can load the .pkl files that pyE2Dpolar creates
aero.bl = load_dict(os.path.join(aero.run_dir, 'bl_dat_prof'))
except:
# run extraction if the BL data is not available
aero.ext_bls(rerun_bl_ext=True, outputname='bl_dat_prof')
# This will also plot some info about the BL in the AoA folder
# as 'bl_data_gl<grid-level>.pdf'. The dict holding the BL data
# stores it for each AoA for which you want to extract the BL
# and saves the data in the run directory as 'bl_dat_prof.pkl'.
# Have a look at the method description to see what is being
# stored, this goes for all stuff really that is being saved.
# Always have a look at the method docstring.
# we can interpolate BL profiles from those extracted using another tool
blq_profiles = interp_bl_prof(xq, aero.bl['gl1']['a2.00']['bl_data'],
aero.bl['gl1']['a2.00']['bl_profiles'])
# plot the interpolated BL profiles just for fun
plt.figure()
viridis = plt.cm.viridis(np.linspace(0, 0.9, len(xq)))
for i in range(len(xq)):
plt.plot(blq_profiles['vt'][i, :] / blq_profiles['U_edge'][i],
blq_profiles['yn'][i, :] / blq_profiles['delta'][i],
# label='$x/c=${:.2f}'.format(blq_profiles['x'][i]),
label='$x/c=${:.2f}'.format(xq[i]),
color=viridis[i], alpha=0.9, lw=1.5)
plt.xlabel(r'$u/U_e$')
plt.ylabel(r'$y/\delta$')
plt.legend()
plt.tight_layout()
plt.savefig('my_bl.pdf')
plt.clf()
# ----- Noise -----
# One reason to extract boundary layers is to compute trailing edge
# noise. Frba and asfi have developed some methods that use the
# BL information extracted from the CFD simulations at the TE to
# determine the noise emissions. This can also be done here. This is
# fast, it is the BL extraction that is slow. The BL is not extracted
# again if not wished for. aero.aoas is used here as well and all noise
# data is stored in the run directory in 'bl_noise.pkl' which then
# includes both noise and BL data.
aero.ext_noise(rerun_bl_ext=False, outputname='bl_noise')
# ======== QUICK METHODS =====================================================
# There are some methods that combine the different elements described above.
# If you want to launch the computations and extract the polars once it is
# finished just use:
# aero.run_extract_polar()
# Or if you do not need any of the inbetween steps just use:
# aero.run_all()
# This will create the mesh, run the simulations and extract the data.
advanced_polars.py
This example script advanced_polars.py is located in examples/. This code presents a possible structure to running several different aerofoils, flow cases and Reynolds numbers and how to use PGL to manipulate the aerofoil surfaces. This is meant as a starting point, however there could many other ways to setup your own polar computations of course.
Note that this code has several checks build in to ensure that not everything (mesh generation etc) is not rerun each time you execute the code. It also launches the computations and does not wait for the the simulations to finish as they can take quite some time to finish. Once the simulations are finished, just relaunch the code and it will extract the polars.
# =====================
# advanced_polars.py
# =====================
# This script is shows how pye2dpolar can easily be adapted to run
# polars for different airfoils and even create erosion like leading
# edge perturbations using PGL. The script should be run twice,
# once for launching all the computations and another time to extract
# all profiles. One could also include other postprocessing tools
# in the loop if wanted, up to you. Be careful with the spectral
# patch generation, it is not fast!
# =====================
import numpy as np
import os
from copy import deepcopy
from e2dpolar.main import e2dpolar
import time
from PGL.components.airfoil import AirfoilShape
import matplotlib.pyplot as plt
from e2dpolar.in_out import load_dict
# ----Run flags
# generate grid from scratch
remake_grid = False
# rerun polar computations
rerun = False
# submit polars to the queue
run = True
# wait in-between launching a number of simulations
wait_finish = 0.1 * 60.
# ----Erosion level and inputs
# note that at level <= 1 no spectral roughness surface is generated
# so there is no dependence on these inputs and you should not iterate
# through them
level = 4
# random seed for phase shifts
seeds = [0]
# ratio of lateral erosion patch to extract (0 <-> 1)
patch_slices = [0.0]
# ----Polar setup
cwd = os.getcwd()
airfoil_files = os.listdir('airfoils')
grid_ni = 512
grid_nj = 256
bsize = 128
# small first grid cell needed for roughness model, which requires y+ < 0.2
grid_first_cell_size = 1.e-7
grid_normal_distribution = [3,
0.0, grid_first_cell_size, 0.0,
1.875e-2, -1, 0.5,
1.0, 2.4, 1.0]
# PLCT turbulence intensity
ti = 0.001
aoas = np.arange(0., 12., 4.)
cpu_time = '00:10:00'
queue = 'workq'
# Reynolds numbers to perform simulations for
# Res = [5e6, 15e6]
Res = [5e6]
# Only perform rough simulations, as surface roughness is applied
# over the eroded region
cases = ['tu']
# ----Runs
# create a folder holding all results
main_folder = 'level' + str(level)
os.makedirs(main_folder, exist_ok=True)
# move into main_folder/
os.chdir(main_folder)
# loop over all airfoil files
for airfoil_file in airfoil_files:
# airfoil name
airfoil = airfoil_file.split('.')[0]
# make airfoil dir
os.makedirs(airfoil, exist_ok=True)
# into airfoil/
os.chdir(airfoil)
print('### running: ' + airfoil + ' ###')
# in this example we are two erosion patterns, that is
# why loop over them here. This is not needed if you
# are running clean aerofoils of course.
for seed, patch_slice in zip(seeds, patch_slices):
# create the run_folders for each seed and slice
runfolder = 'seed{}_nslice{:.2f}'.format(seed, patch_slice)
print(runfolder)
os.makedirs(runfolder, exist_ok=True)
# move into run folder
os.chdir(runfolder)
# ==== CFD sims ====
aero = e2dpolar()
# read airfoil data
aero.pf_coords = np.loadtxt(cwd + '/airfoils/' + airfoil_file)
# sim definitions
aero.grid_ni = grid_ni
aero.grid_nj = grid_nj
aero.grid_bsize = bsize
aero.grid_first_cell_size = grid_first_cell_size
aero.grid_normal_distribution = grid_normal_distribution
# ==== airfoil manipulations ====
# here we use PGL to manipulate the surface. Have a look at
# build/PGL/PGL/components/airfoil.py for all the different options how to
# manipulate airfoil shapes. You can add steps, stallstrips and many
# other nasty things. Here we add some erosion pattern and if you
# go to level 4 also a forward facing step.
aero.af = AirfoilShape(points=aero.pf_coords)
# keep input profile coordinates
aero.af.orig_points = aero.af.points.copy()
# make surface mesh
# this routine also makes sure the coordinates are running in the correct direction,
# normalises them by chord and opens the trailing edge if needed. The opening is
# needed for PGL to create a main_seg with certain attributes need in the erosion
# generation
aero.pf_open_te = 3.e-3
aero.make_surface_mesh(rerun=True)
# save clean airfoil coordinates
pclean = aero.af.points.copy()
if level > 0:
# create eroded LE till x/c < 0.024 or 2.4%
xc = 0.024
rough_patches = {'start': [-xc, 'x'], 'end': [xc, 'x'],
'height': [100e-6, 'ks']}
rough_patches_s = aero._convert_rough_input(rough_patches)
# apply leading edge roughness in eroded region and a little around
aero.rough = True
aero.rough_patches = [deepcopy(rough_patches_s)]
aero.rough_patches[0]['start'][0] -= 20e-3
aero.rough_patches[0]['end'][0] += 20e-3
# -----grid generation
if (not os.path.isfile(os.path.join('grid', 'grid.X2D'))) or remake_grid:
if level > 1:
step = False
if level == 2:
# P2 roughness parameters
stdev = 0.0838e-3
lambda0 = 1.7387e-3
mu = -0.12763e-3
elif level > 2:
# P3 roughness parameters
stdev = 0.2023e-3
lambda0 = 2.3023e-3
mu = -0.4693e-3
if level == 4:
step = True
# From Anders calculations
A = 0.175
B = 0.290
Lx = 203.7447e-3
domega = 2 * np.pi / Lx
# one can store the patch if one likes to do so
patch_name = 'sfgrid/spectralLER_' + airfoil + '_xc{:.3f}_level{}_seed{}'.format(xc, level, seed)
store_patch = True
patch = None
if os.path.isfile(patch_name):
store_patch = False
patch = load_dict(patch_name)
# generate erosion patch and apply it to the aerofoil (check PGL documentation for details)
aero.af.spectralLER(rough_patches_s['start'][0], rough_patches_s['end'][0], ni=256, ds_fac=2.,
patch=patch, store_patch=store_patch, patch_outfilename=patch_name,
domega=domega, fac_dxy_domega=2,
Lx_in=Lx * 0.05,
stdev=stdev, lambda0=lambda0, mu=mu,
damage_limits=[-1.5e-3, 0.],
Nomega=101, Ndirs=21, A=A, B=B, seed=seed,
avgfit_c0=0.4021, avgfit_c1=1.5361,
stdfit_c0=0.3968, stdfit_c1=1.4084, edge_smooth=0.1,
patch_slice=patch_slice,
step=step, s_len_step=40e-3,
h_step=-1.5e-3, h_ni_min=5)
# total number of points after erosion addition
aero.grid_ni = aero.af.ni
# plot some LE roughness output
fig, axs = plt.subplots(1, 1, figsize=(8, 10))
axs.plot(pclean[:, 0], pclean[:, 1], lw=0.5)
axs.plot(aero.af.points[:, 0],
aero.af.points[:, 1], '.-', ms=1, lw=0.5)
axs.set_xlim([-0.01, 0.04])
axs.set_ylim([-0.04, 0.04])
axs.set_aspect('equal', 'box')
fig.savefig('sfgrid/le.pdf')
fig.clf()
aero.make_volume_mesh(rerun=True)
else:
aero.make_mesh(rerun=True)
# -----run polars and extract
if run:
for Re in Res:
for case in cases:
aero.Re = Re
aero.case = case
aero.ti = ti
aero.aoas = aoas
aero.cpu_time = cpu_time
aero.queue = queue
aero.start_glevel = 2
# get the folder name
aero.update()
# submit CFD simulations
launched = aero.run_polar(rerun=rerun)
# when setting rerun=False here, only the aoas are launched
# where no output has been generated, setting it to true all
# earlier results are deleted and rerun
# files to search for which determines if the polars are
# to be extracted
res_file2 = os.path.join(aero.run_dir, 'polar.pkl')
# flag checking whether any CFD simulations were launched
if (not os.path.isfile(res_file2)) and (not launched):
aero.ext_polar()
if launched:
# wait a little before submitting the next one
time.sleep(wait_finish)
os.chdir('../')
# back to main_folder/
os.chdir('../')
Other examples
step.py
Another more basic example how to simulate an aerofoil with a step.
import matplotlib.pyplot as plt
from PGL.components.airfoil import AirfoilShape
from e2dpolar.main import e2dpolar
import numpy as np
# Load instance
aero = e2dpolar()
aero.pf_coords = np.loadtxt('airfoils/naca63-418.dat')
aero.grid_ni = 512
aero.grid_nj = 256
aero.grid_bsize = 128
# ===== SURFACE GRID
# Use the AirfoilShape class from PGL to create a step on the surface of the
# aerofoil. Have a look at PGL to understand how you can manipulate surfaces.
# It is pretty powerful!
# create surface mesh without step
aero.make_surface_mesh(rerun=True)
# store the clean airfoil (ie without step)
pclean = aero.af.points.copy()
# find start and end location in terms of curve length (3.6% x/c according to report)
rough_patches = {'start': [-0.036, 'x'], 'end': [0.036, 'x'],
'height': [1e-4, 'd0']}
rough_patches_s = aero._convert_rough_input(rough_patches)
# step height (positive out of surface)
h = 1e-3
# grid size over step
ds_h = 0.1e-3 / 3
# extra points resolving step, ie added to grid_ni to give the total
ni_step = 128
aero.af.step(rough_patches_s['start'][0], rough_patches_s['end'][0], h, ds_h, ni_step)
# update grid size (original grid + the added points resolving the step)
aero.grid_ni = aero.af.ni
# plot the leading edge with the step on it (modify to make it better)
plt.figure(figsize=[10, 10])
plt.plot(pclean[:, 0], pclean[:, 1], lw=0.5)
plt.plot(aero.af.points[:, 0], aero.af.points[:, 1], '.-', ms=1, lw=0.5)
axs.set_xlim([-0.01, 0.04])
axs.set_ylim([-0.04, 0.04])
axs.set_aspect('equal', 'box')
plt.savefig('sfgrid/le.pdf')
aero.Re = 3e6
# aero.case = 'tr' # transitional boundary layer using e^N transition model
aero.case = 'tu' # fully turbulent boundary layer
# add some roughness over the elevated region to simulate sandpaper
aero.rough = True # switch on the roughness model
aero.rough_patches = [rough_patches_s]
aero.rough_else = {'height': [1.e-25, 'd0']} # roughness outside rough region (needed otherwise the model fails)
aero.ti = 0.001
aero.aoas = np.arange(-4., 4., 2.)
aero.cpu_time = '00:10:00'
aero.procs = 8 # no. of processors (default: 1), if you like it fast
aero.queue = 'workq'
aero.start_glevel = 1
aero.make_volume_mesh(rerun=True)
aero.run_polar(rerun=False, dry_run=False)
aero.check_finish()
aero.ext_polar()
ext_fields = ['u', 'v', 'p']
aero.ext_flowfield(2.0, ext_fields=ext_fields, plt_fields=ext_fields,
write_plot3d=False)
Inspecting grids
Its never a bad idea to check whether your computational grid is actually as you desired it to be. For that purpose some small tools are coming with pye2dpolar that should help you with that, they are located in tools/ but on installation an alias was written to your .bashrc that simplifies it. Just move to the grid/ folder inside the run directory which contains the grid files and type:
pltgrid
which will plot the grid in hpglview. Note that you need to be on the login node for the image to display (unless you have enables -X forwarding from the interactive node). You can zoom in by holding the right mouse button and drawing a zoom box and then using the mid mouse button (three fingers for touchpads) to zoom in. Double click the middle button to zoom back, exit by holding ctrl+c. It is primitive but does the job.
Paraview
Kenneth Lønbaek developed a EllipSys plugin for Paraview with which you can directly read in
*.RST.0* restart files, no need to run the postprocessor. Install the latest version of
Paraview and get the plugin: https://gitlab.windenergy.dtu.dk/kenloen/ellipsys_paraview_plugin.
Follow the documentation for installation and advanced ways of using it.
Loading .pkl files
All data is stored in python dictionaries and saved as pickle. To load the data use the following code snippet
import pickle
def load_dict(name):
with open(name + '.pkl', 'rb') as f:
return pickle.load(f)
So you can load the data in python using:
dict = load_dict(filename)
If you have pye2dpolar installed in your current environment, you could also just import the method:
from e2dpolar.in_out import load_dict
polar = load_dict('polar')