POD from pointclouds

We will now proceed to explain how to perform POD from point clouds. In this instance, we test only for POD in serial, as to perform in parallel, a parallel reader/writer is needed.

If you have saved information in hdf5 and have habilitated mpi4py compilation of it, then you could use this code in parallel.

Please note that to run this example, you need to have a set of files. For this one, we use the files produced from example 4-interpolation/6-interpolating_file_sequences.ipynb

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
import h5py

# 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'

Set up the input parameters

[2]:

file_sequence = [f"../4-interpolation/interpolated_fields{str(1+i).zfill(5)}.hdf5" for i in range(0, 5)]
pod_fields = ["u", "v", "w"]
mesh_fname = "../4-interpolation/points.hdf5"
mass_matrix_fname = "../4-interpolation/points.hdf5"
mass_matrix_key = "mass"
k = len(file_sequence)
p = len(file_sequence)

Call the pysemtools routines

[3]:

from pysemtools.rom.phy_pod_wrappers import pod_from_files
# Output
from pyevtk.hl import gridToVTK

pod, ioh, field3d_shape = pod_from_files(comm, file_sequence, pod_fields, mass_matrix_fname, mass_matrix_key, k, p)

Write out the data

In this case we can easily write the modes from vtk

[4]:

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

from pysemtools.rom.phy_pod_wrappers import write_3dfield_to_file

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

# Write out 5 modes
write_3dfield_to_file("phys_pod.vtk", x, y, z, pod, ioh, modes=[i for i in range(0, 5)], field_shape=field3d_shape, field_names=pod_fields)

# Reconstruct all the snapshots with 5 modes
write_3dfield_to_file("phys_pod.vtk", x, y, z, pod, ioh, modes=[i for i in range(0, 5)], field_shape=field3d_shape, field_names=pod_fields, snapshots=[i for i in range(0, len(file_sequence))])

Writing ./phys_pod_mode0.vtk
Writing ./phys_pod_mode1.vtk
Writing ./phys_pod_mode2.vtk
Writing ./phys_pod_mode3.vtk
Writing ./phys_pod_mode4.vtk
Writing ./phys_pod_reconstructed_data_0
Writing ./phys_pod_reconstructed_data_1
Writing ./phys_pod_reconstructed_data_2
Writing ./phys_pod_reconstructed_data_3
Writing ./phys_pod_reconstructed_data_4

Save the data as numpy arrays

[5]:

from pysemtools.rom.phy_pod_wrappers import save_pod_state

save_pod_state("pod_state.hdf5", pod)

Perform tests

[6]:

from pysemtools.rom.phy_pod_wrappers import physical_space

# Load the mass matrix
with h5py.File(mass_matrix_fname, 'r') as f:
    bm_v = f["mass"][:]

# Go through the snapshots in the file sequence
total_snapshot_energy = 0
for j, fname in enumerate(file_sequence):
    print(f"Testing snapshot {j}: {fname}")

    with h5py.File(fname, 'r') as f:

        # Check one snapshot at a time
        # Here one could also just use the phys that has been previously computed.
        phys = physical_space(pod, ioh, modes=[i for i in range(0, k)], field_shape=field3d_shape, field_names=pod_fields, snapshots=[j])

        # Check if the reconstruction is accurate
        for field in pod_fields:
            passed = np.allclose(phys[j][field], f[field][:])
            print(f"Field {field} passed: {passed}")

        # Check if the energy was accurately captured
        snapshot_energy = 0
        for field in pod_fields:
            snapshot_energy += np.sum(f[field][:]**2*bm_v)
        total_snapshot_energy += snapshot_energy

        reconstruction_energy = 0
        a_i = np.diag(pod.d_1t)@pod.vt_1t[:,j]
        reconstruction_energy += np.sum(np.abs(a_i)**2)

        print(f"Snapshot energy: {snapshot_energy}")
        print(f"Reconstruction energy: {reconstruction_energy}")

        print("=======================================")

# Check if the total energy is kept in the singular values
print(f"total snapshot energy: {total_snapshot_energy}")

total_pod_energy = 0
total_pod_energy += np.sum(pod.d_1t**2)
print(f"Total POD energy: {total_pod_energy}")
print("=======================================")

Testing snapshot 0: ../4-interpolation/interpolated_fields00001.hdf5
Field u passed: True
Field v passed: True
Field w passed: True
Snapshot energy: 3.698958237844187e-05
Reconstruction energy: 3.698959791995535e-05
=======================================
Testing snapshot 1: ../4-interpolation/interpolated_fields00002.hdf5
Field u passed: True
Field v passed: True
Field w passed: True
Snapshot energy: 3.889881492496249e-05
Reconstruction energy: 3.8898829611720744e-05
=======================================
Testing snapshot 2: ../4-interpolation/interpolated_fields00003.hdf5
Field u passed: True
Field v passed: True
Field w passed: True
Snapshot energy: 4.153454569587415e-05
Reconstruction energy: 4.1534561226701525e-05
=======================================
Testing snapshot 3: ../4-interpolation/interpolated_fields00004.hdf5
Field u passed: True
Field v passed: True
Field w passed: True
Snapshot energy: 4.442507641541247e-05
Reconstruction energy: 4.442509363710413e-05
=======================================
Testing snapshot 4: ../4-interpolation/interpolated_fields00005.hdf5
Field u passed: True
Field v passed: True
Field w passed: True
Snapshot energy: 4.68831721967345e-05
Reconstruction energy: 4.68831941768144e-05
=======================================
total snapshot energy: 0.0002087311916114255
Total POD energy: 0.00020873127657229616
=======================================