Generating structured meshes

In most cases, POD can be done directly on a SEM mesh. This of course allows to use derivatice routines present in this framwework to further produce reduced order models.

In some particular cases, however, having an structured mesh might be beneficial. For this, we give some tools in this notebook.

This example might have some overlapping topics with example 4.1

Import general modules

[1]:

# 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

# This example is designed to work in one rank only
if comm.Get_size() > 1:
    raise ValueError("This example is designed to run with one rank only")

Import modules from pysemtools

[2]:

from pysemtools.io.ppymech.neksuite import pynekread, pynekwrite

Create an structured mesh

It is as simple as writing using numpy meshgrid functionality

In this case we want to get a cylinder

[3]:

import pysemtools.interpolation.pointclouds as pcs

# Generate the bounding box of the points
r_bbox = [0, 0.05]
th_bbox = [0, 2*np.pi]
z_bbox = [0 , 1]
nx = 30
ny = 30
nz = 80

# Generate the 1D mesh
r_1d = pcs.generate_1d_arrays(r_bbox, nx, mode="equal")
th_1d = pcs.generate_1d_arrays(th_bbox, ny, mode="equal")
z_1d = pcs.generate_1d_arrays(z_bbox, nz, mode="equal", gain=1)

# Generate differetials (dr, dth, dz)
dr_1d  = pcs.generate_1d_diff(r_1d)
dth_1d = pcs.generate_1d_diff(th_1d, periodic=True) # This is needed to give the same weight to the first and last points as for the other ones. Needed if fourier transform will be applied.
dz_1d  = pcs.generate_1d_diff(z_1d)

# Generate a 3D mesh
r, th, z = np.meshgrid(r_1d, th_1d, z_1d, indexing='ij')
# Generate 3D differentials
dr, dth, dz = np.meshgrid(dr_1d, dth_1d, dz_1d, indexing='ij')

# Generate xy coordinates, which are needed for probes
x = r*np.cos(th)
y = r*np.sin(th)

Mass matrix

To perform POD, we will need a mass matrix. For this we must find a way to combine the coordinates to get the volume of our domain

In our case, because we have cylinder, we now that the volume is calculated by:

\[V = \int_{z_0}^{z_1} \int_{\theta_0}^{\theta_1} \int_{r_0}^{r_1} r \, dr \, d\theta \, dz\]

Therefore, our mass matrix will simply be the terms that are multiplied in the integral:

\[B = r \, dr \, d\theta \, dz\]

If one wants the area for a given angle, then the integral to use is.

\[A = \int_{z_0}^{z_1} \int_{r_0}^{r_1} dr \, dz\]

[ ]:

# For volume
B = r * dr * dth * dz
# For area given each angle slice
A = dr * dz

# Verify that the volume of the cylinder is correct
print(np.sum(B))
print(np.sum(A[:,0,:]))
print(r_bbox[1]**2*np.pi*z_bbox[1])
print(r_bbox[1]*z_bbox[1])

Another alternative is to of course interpolate the mass matrix of a SEM mesh into this set of points, if you have it. But we leave that for another time

Writing the data

If one wishes to use this information in the probes module, then the xyz coordinates as a list will be needed. So one can create that as well:

[5]:

import pysemtools.interpolation.utils as interp_utils

xyz = interp_utils.transform_from_array_to_list(nx,ny,nz,[x, y, z])

Writing as csv

One option is to write everything as a CSV file. In this instance, the data needs to be arrayed as a list already.

[6]:

fname = 'points.csv'
with open(fname, 'w') as f:
    for i in range((xyz.shape[0])):
        f.write(f"{xyz[i][0]},{xyz[i][1]},{xyz[i][2]}\n")

Writing in HDF5

Another flexible option is to write an HDF5 file. In this case one can include an array list xyz, but also more info.

[7]:

import h5py

fname = 'points.hdf5'
with h5py.File(fname, 'w') as f:

    # Create a header
    f.attrs['nx'] = nx
    f.attrs['ny'] = ny
    f.attrs['nz'] = nz
    #f.attrs['probe_list_key'] = 'xyz'

    # Include data sets
    f.create_dataset('x', data=x)
    f.create_dataset('y', data=y)
    f.create_dataset('z', data=z)
    f.create_dataset('r', data=r)
    f.create_dataset('th', data=th)
    f.create_dataset('mass', data=B)
    f.create_dataset('mass_area', data=A)
    #f.create_dataset('xyz', data=xyz)

the interpolation functions will try to find the keyword probe_list_key and use the associated value as probes.

If the probe_key_list does not exist, it will try to use the key xyz by default. If this one also does not exist, it will assemble the list form the x, y, z values

[ ]:

print(np.allclose(B[:,0,:], B[:,1,:]))
print(np.allclose(B[:,2,:], B[:,1,:]))
print(np.allclose(B[:,0,:], B[:,-1,:]))
print(np.sum(B[:,0,:])*ny)