Usage¶

Pymech is a simple interface to read / to write Nek5000 and SIMSON specific data files to / from the Python world. With this capability you could:

make publication-quality figures
post-process
manipulate meshes
generate initial conditions
interpolate or extrapolate solution fields, and

possibly many more potential use-cases – the limit is your imagination! Here we look at some simple operations you can do. Start by installing pymech and:

!git clone --recursive https://github.com/eX-Mech/pymech.git
%cd pymech/tests/data/nek

Cloning into 'pymech'...
remote: Enumerating objects: 3166, done.
remote: Counting objects: 100% (807/807), done.
remote: Compressing objects: 100% (288/288), done.
Receiving objects:   7% (222/3166)
Receiving objects:  25% (792/3166)
remote: Total 3166 (delta 625), reused 557 (delta 512), pack-reused 2359 (from 1)
Receiving objects: 100% (3166/3166), 12.88 MiB | 36.84 MiB/s, done.
Resolving deltas:  21% (449/2135)
Resolving deltas: 100% (2135/2135), done.
Submodule 'data' (https://github.com/eX-Mech/pymech-test-data.git) registered for path 'tests/data'
Cloning into '/tmp/tmpv1mljgpf/pymech/tests/data'...
remote: Enumerating objects: 240, done.        
remote: Counting objects: 100% (240/240), done.        
remote: Compressing objects:  10% (16/157)
remote: Compressing objects:  11% (18/157)
remote: Compressing objects:  12% (19/157)
remote: Compressing objects:  13% (21/157)
remote: Compressing objects:  14% (22/157)
remote: Compressing objects:  15% (24/157)
remote: Compressing objects:  16% (26/157)
remote: Compressing objects:  18% (29/157)
remote: Compressing objects:  19% (30/157)
remote: Compressing objects:  21% (33/157)
remote: Compressing objects: 100% (157/157), done.        
Receiving objects:  31% (75/240)
Receiving objects:  32% (77/240)
Receiving objects:  34% (82/240)
Receiving objects:  35% (84/240)
Receiving objects:  36% (87/240), 16.41 MiB | 32.80 MiB/s
Receiving objects:  37% (89/240), 16.41 MiB | 32.80 MiB/s
remote: Total 240 (delta 82), reused 239 (delta 81), pack-reused 0 (from 0)        
Receiving objects:  45% (108/240), 16.41 MiB | 32.80 MiB/s
Receiving objects: 100% (240/240), 33.88 MiB | 37.10 MiB/s, done.
Resolving deltas: 100% (82/82), done.
Submodule path 'tests/data': checked out '9908c266c7d28ee667a6337ee3cfe73ba577fe3e'
/tmp/tmpv1mljgpf/pymech/tests/data/nek

/home/docs/checkouts/readthedocs.org/user_builds/pymech/envs/stable/lib/python3.9/site-packages/IPython/core/magics/osm.py:417: UserWarning: using dhist requires you to install the `pickleshare` library.
  self.shell.db['dhist'] = compress_dhist(dhist)[-100:]

`pymech.neksuite`¶

from pymech.neksuite import readnek

field = readnek('channel3D_0.f00001')

Simply typing the read field would give you some basic information:

field

<pymech.core.HexaData>
Dimensions:    3
Precision:     4 bytes
Mesh limits:
  * x:         (0.0, 6.2831854820251465)
  * y:         (-1.0, 1.0)
  * z:         (0.0, 3.1415927410125732)
Time:
  * time:      0.2
  * istep:     10
Elements:
  * nel:       512
  * elem:      [<elem centered at [ 0.39269908 -0.98        0.19634954]>
                ...
                <elem centered at [5.89048618 0.98       2.94524309]>]

`pymech.core.HexaData`¶

The pymech.core.HexaData class is the in-memory data structure widely used in Pymech. It stores a description of the field and other metadata. Let’s look at the available attributes:

[attr for attr in dir(field) if not attr.startswith('__')]

['check_bcs_present',
 'check_connectivity',
 'elem',
 'elmap',
 'endian',
 'get_points',
 'istep',
 'lims',
 'lr1',
 'merge',
 'nbc',
 'ncurv',
 'ndim',
 'nel',
 'offset_connectivity',
 'time',
 'update_ncurv',
 'var',
 'wdsz']

Here check_connectivity and merge are methods, and elem and elmap are data objects. The rest of the attributes store metadata. See pymech.core to understand what they imply.

print(field.endian, field.istep, field.lr1, field.nbc, field.ncurv, field.ndim, field.nel, field.time, field.var, field.wdsz)

little 10 (8, 8, 8) 0 [] 3 512 0.2 (3, 3, 1, 0, 0) 4

The elem attribute contains data of physical fields. It is an array of lists, with each array representing an element.

print("There are", field.nel, "elements in this file")

There are 512 elements in this file

`pymech.core.Elem`¶

The raw arrays are stored as a list of pymech.core.Elem instances as pymech.core.HexaData’s elem attribute.

first_element = field.elem[0]
print("Type =", type(first_element))
print(first_element)

Type = <class 'pymech.core.Elem'>
<elem centered at [ 0.39269908 -0.98        0.19634954]>

Let us look at the attributes of an element

[attr for attr in dir(first_element) if not attr.startswith('__')]

['bcs',
 'ccurv',
 'centroid',
 'curv',
 'face_center',
 'pos',
 'pres',
 'scal',
 'smallest_edge',
 'temp',
 'vel']

Except for the following attributes of an element object

print(first_element.bcs, first_element.ccurv, first_element.curv)

[] ['', '', '', '', '', '', '', '', '', '', '', ''] [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

it contains large arrays

print("Shape of element velocity and pressure arrays = ", first_element.vel.shape, first_element.pres.shape)

Shape of element velocity and pressure arrays =  (3, 8, 8, 8) (1, 8, 8, 8)

`pymech.dataset`¶

from pymech.dataset import open_dataset

ds = open_dataset('channel3D_0.f00001')

This function loads the field file in a more convenient xarray dataset.

The dataset is more descriptive and useful for exploratory work, such as post-processing and plotting.

Visualizing¶

Average in spanwise (z) direction and plot velocity profile

ds.ux.mean('z').plot()

<matplotlib.collections.QuadMesh at 0x7f713df21550>

_images/8a6a2eb6473b5259a3170d8312e43333218e9ead79d8739723c7225db9aff6ca.png

Average in both horizontal direction and plot 1D profile

ds_1d = ds.mean(['x', 'z'])
ds_1d.ux.plot()

[<matplotlib.lines.Line2D at 0x7f71400caeb0>]

_images/d731c764e0449ce55761526fd66c1c93ae67e41ee5b972fe5a2aa01101339285.png

It is also worth knowing that it is possible to:

Parallelize these operations using ds.chunk method followed by ds.compute
Open a multiple files into a single dataset using pymech.dataset.open_mfdataset, optionally in parallel.

Read the xarray documentation to see how to use them.

`pymech.meshtools`¶

This class helps working with Nek5000 meshes. It can modify existing meshes and generate new ones.

2D Mesh generation¶

Boxes¶

You can generate a rectangular box mesh like this:

import pymech.meshtools as mt

# box of size 5×1
xmin = 0
xmax = 5
ymin = 0
ymax = 1
# 20 × 10 elements resolution
nx = 20
ny = 10
# velocity inflow on the left
bc_inflow = ['v']
# outflow on the right
bc_outflow = ['O']
# wall at the bottom
bc_wall = ['W']
# free-stream (fixed vx and d(vy)/dy = 0) at the top
bc_freestream = ['ON']
box = mt.gen_box(
    nx,
    ny,
    xmin,
    xmax,
    ymin,
    ymax,
    bcs_xmin=bc_inflow,
    bcs_xmax=bc_outflow,
    bcs_ymin=bc_wall,
    bcs_ymax=bc_freestream,
)

It is also possible to start from a square box and map it to a different shape, or simply move the elements in the box to change the resolution locally. For example, this refines the resolution close to the wall such that the first element has a height of 0.025 instead of 0.1.

l0 = 0.025
alpha = mt.exponential_refinement_parameter(l0, ymax, ny)

# maps (x, y) -> (x', y')
def refinement_function(x, y):
    iy = ny * y / ymax
    return (x, l0 * (1 - alpha**iy)/(1 - alpha))

# This mapping doesn't introduce any curvature, so we don't need to encode any new curvature in the mesh
box = mt.map2D(box, refinement_function, curvature=False, boundary_curvature=False)

Other geometries¶

Pymech also lets you generate a disc mesh. The following example generates a disc domain with a temperature field

# the disc will have radius 1
radius = 1
# centre box of 5×5 elements
n_centre = 5
# four side boxes of 6×5 elements
n_bl = 6
# the diagonal of the centre box will be half of the disc diametre 
s_param = 0.5
# (x, y), (vx, vy), p, t, no passive scalars
variables = [2, 2, 1, 1, 0]
# boundary conditions: adiabatic wall
bcs_sides = ['W', 'I']
# generate the mesh
disc_mesh = mt.gen_circle(radius, s_param, n_centre, n_bl, var=variables, bc=bcs_sides)

Extrusion¶

Pymech can extrude 2D meshes into 3D ones.

We can for example extrude the disc into a cylinder with velocity and temperature inlet at the bottom and outflow at the top:

# extrude in ten elements vertically between -1 and +1
z = [-1, -0.9, -0.7, -0.5, -0.25, 0, 0.25, 0.5, 0.7, 0.9, 1]
bc_inflow = ['v', 't']
bc_outflow = ['O', 'I']
# bc1 denotes the boundary conditions at z = -1, and bc2 at z = +1
cylinder_mesh = mt.extrude(disc_mesh, z, bc1=bc_inflow, bc2=bc_outflow)

Usage¶

`pymech.neksuite`¶

`pymech.core.HexaData`¶

`pymech.core.Elem`¶

`pymech.dataset`¶

Computing statistics¶

Slicing¶

Visualizing¶

`pymech.meshtools`¶

2D Mesh generation¶

Boxes¶

Other geometries¶

Extrusion¶

Usage¶

pymech.neksuite¶

pymech.core.HexaData¶

pymech.core.Elem¶

pymech.dataset¶

Computing statistics¶

Slicing¶

Visualizing¶

pymech.meshtools¶

2D Mesh generation¶

Boxes¶

Other geometries¶

Extrusion¶

`pymech.neksuite`¶

`pymech.core.HexaData`¶

`pymech.core.Elem`¶

`pymech.dataset`¶

`pymech.meshtools`¶