UQ of velocity fields

In this notebook we show how one can use multiple batches of averages fields to perform a UQ analysis.

For this case we assume that you have a set of files already. The visualizations here rely on you having an structured mesh

Import general modules

# Import required modules
from mpi4py import MPI #equivalent to the use of MPI_init() in C
import matplotlib.pyplot as plt
import numpy as np

# Get mpi info
comm = MPI.COMM_WORLD

# Hide the log for the notebook. Not recommended when running in clusters as it is better you see what happens
import os
os.environ["PYSEMTOOLS_HIDE_LOG"] = 'true'

Define inputs

query_points_fname = "./points.hdf5" # In our case this file contains the mesh and the mass matrix
file_sequence = [f"mean_field{str(i+1).zfill(5)}.hdf5" for i in range(0, 10)]

Obtain the mean and variance

For this case we can simply use the Non overlapped batch mean method, assuming that we have a number of batch means.

from pysemtools.io.wrappers import read_data
from pysemtools.postprocessing.statistics.uq import NOBM

# Read the mesh dat for latter plotting
mesh_data = read_data(comm, query_points_fname, ["x", "y", "z"], dtype = np.single)

# Compute the mean and variance of the specified fields
mean, var = NOBM(comm, file_sequence, ["u", "v", "w"])

Plot it

The mean

from pysemtools.postprocessing.plotting import isosurfaces
from IPython.display import Image, display

isos = {}
#isos["u"] = np.linspace(-0.007, 0.007, 2)
#isos["v"] = np.linspace(-0.007, 0.007, 2)
#isos["w"] = np.linspace(-0.015, 0.015, 2)

isos["u"] = np.linspace(np.min(mean["u"]), np.max(mean["u"]), 4)
isos["v"] = np.linspace(np.min(mean["v"]), np.max(mean["v"]), 4)
isos["w"] = np.linspace(np.min(mean["w"]), np.max(mean["w"]), 4)
pl = isosurfaces(mesh_data, mean, isosurfaces = isos, shape = (1, 3), window_size = [1920, 1080])

for plotter in pl:
    # Capture the plot as an image and show it
    image_path = f"static_plot_SME.png"
    plotter.screenshot(image_path)
    plotter.close()
    display(Image(filename=image_path))

The variance

isos = {}
isos["u"] = np.linspace(np.min(var["u"]), np.max(var["u"]), 4)
isos["v"] = np.linspace(np.min(var["v"]), np.max(var["v"]), 4)
isos["w"] = np.linspace(np.min(var["w"]), np.max(var["w"]), 4)
pl = isosurfaces(mesh_data, var, isosurfaces = isos, shape = (1, 3), window_size = [1920, 1080], colormap="cool", opacity=1)

for plotter in pl:
    # Capture the plot as an image and show it
    image_path = f"static_plot_SME.png"
    plotter.screenshot(image_path)
    plotter.close()
    display(Image(filename=image_path))

The interpolated fields and also the fields read from nek5000 have a set of given keys. One could know what each of the keys represent by heart, but we have also provided some “indices” to help identify what is what.

It can be used as folows:

from pysemtools.postprocessing.statistics.file_field_index import file_field_index

# Specify the source of the data
key_index = file_field_index(source_code="Neko_old") # Neko old alludes to statistics written between 2 files
print(key_index["U"]["pynek_field_key"])

u

perform the rotations and write to file

from pysemtools.postprocessing.rotations import cartesian_to_cylindrical_rotation_matrix as rotation_matrix
from pysemtools.postprocessing.rotations import rotate_tensor
import os

# Specify the components of the tensor. This should match the key in the file
for i, fname in enumerate(file_sequence):
    tensor_components = [key_index["U"]["pynek_field_key"], key_index["V"]["pynek_field_key"], key_index["W"]["pynek_field_key"]]
    prefix = os.path.basename(fname).split(".")[0]
    rotate_tensor(comm, query_points_fname, fname, tensor_components, rotation_matrix, output_file=f"rotated_{prefix}.hdf5", output_field_names=["U_r", "U_th", "U_z"])

now apply the NOBM method

rotated_file_sequence = [f"rotated_{os.path.basename(fname).split('.')[0]}.hdf5" for fname in file_sequence]
mean, var = NOBM(comm, rotated_file_sequence, ["U_r", "U_th", "U_z"])

print(np.min(mean["U_r"]), np.max(mean["U_r"]))

-0.011690581479526809 0.04689160678188811

plot the mean

isos = {}
isos["U_r"] = np.linspace(-0.007, 0.007, 2)
isos["U_th"] = np.linspace(-0.007, 0.007, 2)
isos["U_z"] = np.linspace(-0.015, 0.015, 2)
pl = isosurfaces(mesh_data, mean, isosurfaces = isos, shape = (1, 3), window_size = [1920, 1080])

for plotter in pl:
    # Capture the plot as an image and show it
    image_path = f"static_plot_SME.png"
    plotter.screenshot(image_path)
    plotter.close()
    display(Image(filename=image_path))

now the variance

isos = {}
isos["U_r"] = np.linspace(0, 1e-4, 3)
isos["U_th"] = np.linspace(0, 1e-4, 3)
isos["U_z"] = np.linspace(0, 2e-4, 3)
pl = isosurfaces(mesh_data, var, isosurfaces = isos, shape = (1, 3), window_size = [1920, 1080], colormap="cool", opacity=1)

for plotter in pl:
    # Capture the plot as an image and show it
    image_path = f"static_plot_SME.png"
    plotter.screenshot(image_path)
    plotter.close()
    display(Image(filename=image_path))

In this case, we can average in the azimuthal direction.

import h5py

# Load mass matrix
with h5py.File(query_points_fname, 'r') as f:
    bm = f["mass"][:]
bm[np.where(bm == 0)] = 1e-12

# Load the mesh
with h5py.File(query_points_fname, 'r') as f:
    r = f["r"][:]
    th = f["th"][:]
    z = f["z"][:]

u = mean["U_r"]
v = mean["U_th"]
w = mean["U_z"]

#Average in the azimuthal direction direction
u_rz = np.sum(u*bm, axis=1)/ np.sum(bm, axis=1)
v_rz = np.sum(v*bm, axis=1)/ np.sum(bm, axis=1)
w_rz = np.sum(w*bm, axis=1)/ np.sum(bm, axis=1)

levels = np.linspace(-0.05, 0.05, 100)

fig, ax = plt.subplots(1,3, figsize=(20,10))
ax[0].tricontourf(r[:,0,:].flatten(),z[:,0,:].flatten(), u_rz.flatten(), levels=levels, cmap="RdBu_r")
ax[1].tricontourf(r[:,0,:].flatten(),z[:,0,:].flatten(), v_rz.flatten(), levels=levels, cmap="RdBu_r")
ax[2].tricontourf(r[:,0,:].flatten(),z[:,0,:].flatten(), w_rz.flatten(), levels=levels, cmap="RdBu_r")
#plt.show()
fig.savefig("contours.png")
plt.close()
display(Image(filename="contours.png"))

Plot the rotated values in 1d, by averaging in azimuthal and radial directions. In this case, this might not be so telling, as the flow is not homogenous in the radial direction, but it have give an indication on how the flow behaves.

In this case, we might have had too few points close to the wall, as we do not observe any peak in the z velocity component. Note that it would appear better if we take the square root.

# Load mass matrix
with h5py.File(query_points_fname, 'r') as f:
    bm = f["mass"][:]
bm[np.where(bm == 0)] = 1e-12

u = mean["U_r"]
v = mean["U_th"]
w = mean["U_z"]

# Average in the azimuthal direction and radial direction
u_z = np.sum(u*bm, axis=(0,1))/ np.sum(bm, axis=(0,1))
v_z = np.sum(v*bm, axis=(0,1))/ np.sum(bm, axis=(0,1))
w_z = np.sum(w*bm, axis=(0,1))/ np.sum(bm, axis=(0,1))


fig, ax = plt.subplots(1,1, figsize=(20,3))
ax.plot(z[0,0,:], u_z, '-k', label="u_r")
ax.plot(z[0,0,:], v_z, '-b', label="u_th")
ax.plot(z[0,0,:], w_z, '-r', label="u_z")
ax.set_xlabel("z")
fig.legend()
#plt.show()
fig.savefig("mean_fields.png")
plt.close()
display(Image(filename="mean_fields.png"))