{ "cells": [ { "cell_type": "markdown", "id": "5ac68bf1", "metadata": {}, "source": [ "# Shallow Recurrent Decoders (SHRED) for Parametric State Estimation\n", "This notebook demostrates how to combine the *pyforce* package with the [pySHRED](https://github.com/PyShred-Dev/PyShred) and our vanilla implemtation in [NuSHRED](https://github.com/ERMETE-Lab/NuSHRED).\n", "\n", "This ML technique is based on the work by [Williams et al. (2024)](https://royalsocietypublishing.org/rspa/article/480/2298/20240054/66770/Sensing-with-shallow-recurrent-decoder), and the extension to parametric problems by [Tomasetto et al. (2025)](https://www.nature.com/articles/s41467-025-65126-y) and [Riva et al. (2025)](https://arxiv.org/abs/2503.08904). \n", "\n", "Two versions of SHRED are available: the official PySHRED package, and the in-house NuSHRED code. The former can be installed via pip:\n", "\n", "```bash\n", "pip install pyshred\n", "```\n", "The latter can be cloned from the GitHub repository and the code can be used directly:\n", "\n", "```bash\n", "git clone https://github.com/ERMETE-Lab/NuSHRED.git\n", "```\n", "\n", "In this notebook, we will consider a parametric flow over cylinder problem, from CFDBench [Luo et al. (2024)](https://arxiv.org/abs/2310.05963): the inlet boundary is parametrized with respect to its intensity, from 0.1 m/s to 1.0 m/s." ] }, { "cell_type": "markdown", "id": "588b26cf", "metadata": {}, "source": [ "At first, let us load the raw snapshots, downloaded from [the benchmark](https://huggingface.co/datasets/chen-yingfa/CFDBench/blob/main/cylinder/bc.zip) - the zenodo folder only contains the first 10 parameters, instead of the full 50." ] }, { "cell_type": "code", "execution_count": 27, "id": "bddf1cff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of snapshots: 10, Number of time steps: 1000, Grid size: 64 x 64\n" ] } ], "source": [ "import numpy as np\n", "import os\n", "from IPython.display import clear_output\n", "\n", "folder_to_save = '../Datasets/CFDBenchFlowCyl'\n", "\n", "id_cases = [f for f in os.listdir(folder_to_save+'/') if os.path.isdir(os.path.join(folder_to_save+'/', f))]\n", "id_cases.sort()\n", "\n", "snap_data = list()\n", "for case_i in range(1, len(id_cases)):\n", " snap_data.append(dict())\n", "\n", " snap_data[case_i-1]['u'] = np.load(folder_to_save + '/' + id_cases[case_i] + '/u.npy')\n", " snap_data[case_i-1]['v'] = np.load(folder_to_save + '/' + id_cases[case_i] + '/v.npy')\n", "\n", "Nt, Nhx, Nhy = snap_data[0]['u'].shape\n", "Ns = len(snap_data)\n", "\n", "var_names = ['u', 'v']\n", "\n", "print(f'Number of snapshots: {Ns}, Number of time steps: {Nt}, Grid size: {Nhx} x {Nhy}')" ] }, { "cell_type": "markdown", "id": "fa72c541", "metadata": {}, "source": [ "From the benchmark description, the geometry is generated in the pyvista format." ] }, { "cell_type": "code", "execution_count": 28, "id": "1f7c52ea", "metadata": {}, "outputs": [], "source": [ "radius = 0.01\n", "length = 0.16 + 0.06 + 2*radius\n", "height = 2 * 0.06 + 2*radius\n", "\n", "import pyvista as pv\n", "\n", "x = np.linspace(0, length, Nhx)\n", "y = np.linspace(0, height, Nhy)\n", "X, Y = np.meshgrid(x, y, indexing=\"ij\")\n", "points = np.c_[X.ravel(), Y.ravel(), np.zeros_like(X.ravel())]\n", "\n", "# --- Build connectivity for quadrilateral cells ---\n", "# Each quad has 4 vertices, and we store them in the format:\n", "# [4, id0, id1, id2, id3] for each cell\n", "cells = []\n", "for i in range(Nhx - 1):\n", " for j in range(Nhy - 1):\n", " p0 = i * Nhy + j\n", " p1 = p0 + 1\n", " p2 = p0 + Nhy + 1\n", " p3 = p0 + Nhy\n", " cells.append([4, p0, p1, p2, p3])\n", "\n", "cells = np.array(cells, dtype=np.int64).ravel()\n", "\n", "# Cell types: VTK_QUAD = 9\n", "celltypes = np.full((Nhx - 1) * (Nhy - 1), 9, dtype=np.uint8)\n", "\n", "# Plotting the grid\n", "grid = pv.UnstructuredGrid(cells, celltypes, points)\n", "nodes = grid.points[:, :2] # Extract only x and y coordinates" ] }, { "cell_type": "markdown", "id": "2d97d45b", "metadata": {}, "source": [ "The snapshots loaded before are the $x$ and $y$ components of the velocity field, parameterised with respect to the inlet velocity. We are going to store them in numpy arrays of the following shape $N_{params}\\times N_t\\times \\mathcal{N}_h$, where $N_{params}$ is the number of parameters (inlet velocities), $N_t$ is the number of time steps, and $\\mathcal{N}_h$ is the number of spatial degrees of freedom." ] }, { "cell_type": "code", "execution_count": 29, "id": "5ba9388b", "metadata": {}, "outputs": [], "source": [ "snapshots_np = {var_name: np.zeros((Ns, Nt, Nhx * Nhy)) for var_name in var_names}\n", "\n", "for i in range(Ns):\n", " for var_name in var_names:\n", " for tt in range(Nt):\n", " snapshots_np[var_name][i, tt, :] = snap_data[i][var_name][tt].ravel(order='F')\n", "\n", "Nh = Nhx * Nhy\n", "\n", "from matplotlib import patches\n", "def create_circle(ax):\n", " radius = 0.01\n", " center = ((0.06+radius/2), (0.06+radius))\n", " circle = patches.Ellipse(center, 2*radius, 2*radius, edgecolor='black', facecolor='white', linewidth=2)\n", " ax.add_patch(circle)" ] }, { "cell_type": "markdown", "id": "2a2cd7ef", "metadata": {}, "source": [ "## pySHRED implementation\n", "We first demonstrate the pySHRED implementation using POD-based compressive training, either directly from the package or adopted the functionalities of *pyforce*." ] }, { "cell_type": "markdown", "id": "35ef0bc6", "metadata": {}, "source": [ "### Directly from pySHRED package\n", "The parametric version has to be loaded." ] }, { "cell_type": "code", "execution_count": 30, "id": "a71c2086", "metadata": {}, "outputs": [], "source": [ "from pyshred import ParametricDataManager, SHRED, ParametricSHREDEngine\n", "\n", "# Initialize ParametricSHREDDataManager\n", "manager_pod = ParametricDataManager(\n", " lags = 25,\n", " train_size = 0.8,\n", " val_size = 0.1,\n", " test_size = 0.1,\n", " )\n", "\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "markdown", "id": "1e0368f2", "metadata": {}, "source": [ "Let us add the different fields, the component $u$ is the one we want to reconstruct, while $v$ is indirectly reconstructed." ] }, { "cell_type": "code", "execution_count": 31, "id": "ce6e01c9", "metadata": {}, "outputs": [], "source": [ "# Measured Field\n", "manager_pod.add_data(\n", " data = snapshots_np['u'],\n", " id = 'u',\n", " random = 3,\n", " compress = 15\n", ")\n", "\n", "# Unobserved Field\n", "manager_pod.add_data(\n", " data = snapshots_np['v'],\n", " id = 'v',\n", " compress = 15\n", ")" ] }, { "cell_type": "markdown", "id": "c2293689", "metadata": {}, "source": [ "Here's a summary of the measurements locations" ] }, { "cell_type": "code", "execution_count": 32, "id": "dea1f293", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| \n", " | data id | \n", "sensor_number | \n", "type | \n", "loc/traj | \n", "
|---|---|---|---|---|
| 0 | \n", "u | \n", "0 | \n", "stationary (random) | \n", "(3644,) | \n", "
| 1 | \n", "u | \n", "1 | \n", "stationary (random) | \n", "(2751,) | \n", "
| 2 | \n", "u | \n", "2 | \n", "stationary (random) | \n", "(3327,) | \n", "