notebooks/coherent_structures_dim_reduction_exercise.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "basic-adapter",
   "metadata": {},
   "source": [
    "![CC](https://i.creativecommons.org/l/by/4.0/88x31.png)\n",
    "\n",
    "This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).\n",
    "\n",
    "# Analyzing coherent structures in flows displaying transonic buffets\n",
    "\n",
    "This exercise consists of two parts, which are meant to be solved in two individual exercise sessions. In the first part, we analyze coherent structures in the surface pressure coefficient of a NACA0012 airfoil in transonic flow conditions by means of a principal component analysis (PCA). In the second part, we analyze the same data using dynamic mode decomposition (DMD). The specific flow conditions in terms of Reynolds number, Mach number, and angle of attack are $Re=10^7$, $Ma=0.75$, and $\\alpha = 4^\\circ$, respectively. The pressure coefficient is defined as:\n",
    "\n",
    "$$\n",
    "  c_p = \\frac{p}{0.5\\rho_\\infty U_\\infty},\n",
    "$$\n",
    "\n",
    "where $p$ is the total pressure, $U_\\infty$ is the free-stream speed and $\\rho_\\infty$ is the free-stream density. The dominant frequency associated with the shock buffet is $f\\approx 38Hz$.\n",
    "\n",
    "To execute this exercise, create an empty script or Jupyter notebook in the *exercise* folder.\n",
    "\n",
    "## Principal component analysis\n",
    "\n",
    "### Loading and masking the dataset\n",
    "\n",
    "The simulation was conducted using the OpenFOAM solver *rhoCentralFoam*. The [pressure](https://www.openfoam.com/documentation/guides/latest/doc/guide-fos-field-pressure.html) function object combined with [surface sampling](https://www.openfoam.com/documentation/guides/latest/doc/guide-fos-sampling-surfaces.html) was used to compute and write the pressure coefficient. Such data are typically written to the directory *postProcessing/surfaces/time/field.csv*. If you download the latest lecture dataset, you find the *surfaces* folder under *datasets/naca0012_buffet/surfaces/*. To access the data, we use the flowTorch [CSVDataloader](https://flowmodelingcontrol.github.io/flowtorch-docs/1.0/flowtorch.data.html#module-flowtorch.data.csv_dataloader). More specifically, we create a loader object using the `CSVDataloader.from_foam_surface()` class method, which indicates the dataloader how the data are organized:\n",
    "```\n",
    "from flowtorch.data import CSVDataloader\n",
    "\n",
    "data = \"../datasets/naca0012_buffet/surface/\"\n",
    "loader = CSVDataloader.from_foam_surface(data, \"total(p)_coeff_airfoil.raw\")\n",
    "```\n",
    "The loader provides access to several attributes of the dataset:\n",
    "\n",
    "- `loader.vertices` returns the x/y/z coordinates of the sample points located on the surface\n",
    "- `loader.write_times` returns a list of available snapshot times formatted as strings, e.g., \"0.01\"\n",
    "- `loader.load_snapshot(field, time)` loads a *field* at the specified *time*; the *time* must be given as string\n",
    "\n",
    "As a first step, load the vertices, normalize them with the chord length $c=0.60105m$, and visualize the $x$ and $y$ components (first and second colum of *loader.vertices*) as scatter plot.\n",
    "\n",
    "The dataset contains both lower and upper surface (pressure and suction side), but we will analyze only the upper side. Therefore, we create a boolean mask to distinguish between points on the lower and upper surface. flowTorch provides [selection tools](https://flowmodelingcontrol.github.io/flowtorch-docs/1.0/flowtorch.data.html#module-flowtorch.data.selection_tools) to create such masks. Use the *mask_box()* function to select only points on the upper surface. Hint: the airfoil is symmetric with respect to the plane $(x, y=0, z)$.\n",
    "```\n",
    "from flowtorch.data import mask_box\n",
    "mask = mask_box(normalized_vertices, lower=[?, ?, -1.0], upper=[?, ?, 1.0])\n",
    "```\n",
    "To apply the mask to a tensor, we can use the [masked_select()](https://pytorch.org/docs/stable/generated/torch.masked_select.html) function. For example, the commands\n",
    "```\n",
    "x = pt.masked_select(normalized_vertices[:, 0], mask)\n",
    "z = pt.masked_select(normalized_vertices[:, 2], mask)\n",
    "```\n",
    "could be used to sub-select the $x$ and $z$ component of points on the upper surface. The number of selected points is equal to the sum over all elements in *mask*, e.g., `n_points = mask.sum().item()`.\n",
    "\n",
    "Next, we analyze the available write times and select a window for the PCA analysis. To do so,\n",
    "\n",
    "- print the first and the last write time; the elements of `loader.write_times` are sorted\n",
    "- compute the time step between the first two snapshots by converting the strings into floats\n",
    "- create a new list of string-encoded times selecting only the elements of `loader.write_times` with $0.0326s \\le t \\le 0.1015s$\n",
    "\n",
    "Based on the number of selected points and snapshots, we can now allocate an empty data matrix, in which individual snapshots will be organized as column vectors:\n",
    "```\n",
    "data_matrix = pt.zeros((?, ?))\n",
    "```\n",
    "To fill up the data matrix, execute the following steps:\n",
    "\n",
    "- loop over the selected times\n",
    "- load the snapshot at time *t* using loader.load_snapshot(\"f\", t)\n",
    "- apply the mask to the snapshot\n",
    "- save the selected data as the corresponding column of the data matrix\n",
    "\n",
    "After the data matrix is complete, compute the temporal mean and variance of each point in the data matrix and plot both fields using [tricontourf](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.tricontourf.html). The function call might look similar to the following line:\n",
    "```\n",
    "plt.tricontourf(x, z, data_matrix.mean(dim=1))\n",
    "```\n",
    "\n",
    "Finally, subtract the mean from each column. Hint: use the *dim* keyword of the [mean()](https://pytorch.org/docs/stable/generated/torch.mean.html) function and [unsqueeze()](https://pytorch.org/docs/stable/generated/torch.unsqueeze.html) the result to convert the row into a column vector, which can be more easily subtracted from each column of the data matrix.\n",
    "\n",
    "### Computing the singular value decomposition\n",
    "\n",
    "We compute the PCA be means of a singular value decomposition (SVD) of the mean-subtracted data matrix. The [svd()](https://pytorch.org/docs/stable/generated/torch.linalg.svd.html?highlight=svd#torch.linalg.svd) function returns three tensors, namely the left singular vectors, the singular values, and the right singular vectors (in that order).\n",
    "```\n",
    "U, s, VH = pt.linalg.svd(data_matrix, full_matrices=False)\n",
    "# the right singular vectors are usually returned transposed\n",
    "V = VH.T\n",
    "```\n",
    "\n",
    "The singular values tell us about the importance of each mode. The elements of *s* are already organized in descending order. Create the following two plots to visualize the importance of the first 50 modes:\n",
    "\n",
    "- a [bar](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.bar.html) plot of the singular values\n",
    "- a bar plot of the normalized singular values, i.e., `s_rel = [si/s.sum().item() for si in s]`\n",
    "\n",
    "### Visualizing modes and mode coefficients\n",
    "\n",
    "The column vectors of *U* correspond to the co-called principal components or proper orthogonal modes. Visualize the first five modes (the first five columns of *U*) as you visualized mean and variance in the first sub-task. What can you say about the first mode?\n",
    "\n",
    "Convert the selected snapshot times into floating point values and plot the first five columns of *V* over time. To identify the frequency associated with each mode, use the [welch](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html) function as follows:\n",
    "```\n",
    "from scipy.signal import welch\n",
    "\n",
    "# times is a list containing the selected write times as floats\n",
    "dt = times[1] - times[0]\n",
    "for i in range(5):\n",
    "    f, power = welch(V[:, i], 1/dt, nfft=len(times)*2)\n",
    "    plt.plot(f, power, label=f\"mode {i+1}\")\n",
    "    print(f\"Maximum amplitude at f={f[np.argmax(power)]:2.4f}Hz\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlim(0, 300)\n",
    "plt.xlabel(r\"$f$ in $Hz$\")\n",
    "plt.ylabel(\"PSD\")\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "Which mode is most representative of the shock buffet?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "original-dating",
   "metadata": {},
   "source": [
    "## Dynamic mode decomposition\n",
    "\n",
    "In the second part of this exercise, we use the same surface pressure data as before but apply dynamic mode decomposition (DMD) to obtain spatial modes that have only a single frequency associated with them. To get started, assemble the data matrix as in the first part of the exercise but do not extract the temporal mean from the snapshots.\n",
    "\n",
    "### DMD implementation\n",
    "\n",
    "Computing the DMD essentially consists of four steps:\n",
    "\n",
    "1. compute the SVD of the data matrix omitting the last snapshot: $ \\mathbf{X} = \\mathbf{U\\Sigma V}^H $\n",
    "2. compute the reduced linear operator $\\tilde{\\mathbf{A}} = \\mathbf{U}^H\\mathbf{X}^\\prime\\mathbf{V\\Sigma}^{-1}$\n",
    "3. compute the eigen decomposition of the reduced linear operator: $\\tilde{\\mathbf{A}} = \\mathbf{W\\Lambda W}^{-1}$\n",
    "4. reconstruct the eigenvectors (modes) of the full linear operator: $\\mathbf{\\Phi} = \\mathbf{X}^\\prime\\mathbf{V\\Sigma}^{-1}\\mathbf{W}$\n",
    "\n",
    "Apply the DMD algorithm to the pressure surface coefficient data matrix. For the implementation of this and the next tasks, the following hints might be helpful:\n",
    "\n",
    "- [torch.linalg.svd](https://pytorch.org/docs/stable/generated/torch.linalg.svd.html); use `full_matrices=False`\n",
    "- [torch.linalg.eig](https://pytorch.org/docs/stable/generated/torch.linalg.eig.html#torch.linalg.eig)\n",
    "- [torch.diag](https://pytorch.org/docs/stable/generated/torch.diag.html)\n",
    "- to convert a tensor of floats into a complex tensor, use the *type* method, e.g., `matrix.type(pt.cfloat)`\n",
    "- to access real and imaginary parts of a complex tensor, use the properties `matrix.real` and `matrix.imag`\n",
    "- compute the complex conjugate transpose of a matrix with the syntax `matrix.conj().T`\n",
    "- compute matrix multiplications using the @ operator\n",
    "\n",
    "To see if the implementation yields sensibe results, plot the eigenvalues in the complex plane. Most eigenvalues should be located on the unit circle. The following code provides most of the lines to create such a plot.\n",
    "\n",
    "```\n",
    "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n",
    "t = pt.linspace(0, 2 * np.pi, 100)\n",
    "ax.plot(pt.cos(t), pt.sin(t), ls=\"--\", color=\"k\", lw=2)\n",
    "ax.scatter(???, ???, marker=\"x\", lw=1, s=20, zorder=7)\n",
    "ax.set_xlim(-1.3, 1.3)\n",
    "ax.set_ylim(-1.3, 1.3)\n",
    "ax.set_xlabel(r\"$Re(\\lambda)$\")\n",
    "ax.set_ylabel(r\"$Im(\\lambda)$\")\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "### Plotting the spectrum\n",
    "\n",
    "We can now determine the frequency associated with each mode and assign a mode amplitude/importance. The frequency corresponds to the imaginary part of the eigenvalues $\\omega_i$ of the time-continuous linear operator. The eigenvalues of discrete and continuous operator are connected by the relation $\\omega_i = \\mathrm{log}(\\lambda_i)/\\Delta t$. Deviding the imaginary part of $\\omega_i$ by $2\\pi$ yields the frequency of the $i$th mode in $Hz$. Use the initial conditions $ \\mathbf{b} = \\mathbf{\\Phi}^{-1}\\mathbf{x}_0 $ as mode amplitudes (ignore the imaginary part).\n",
    "\n",
    "A common way to visualize a spectrum is a [stem plot](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.stem.html). The eigenvalues come as complex conjugate pairs because our data matrix contains only real values. The following code snippet shows how to plot only the positive frequencies:\n",
    "\n",
    "```\n",
    "amplitude = b.abs()\n",
    "freq = omega.imag / (2.0*np.pi)\n",
    "freq_i_pos = (freq > 0).nonzero().flatten()\n",
    "\n",
    "fig, ax= plt.subplots()\n",
    "ax.stem(freq[freq_i_pos].numpy(), amplitude[freq_i_pos].numpy(), basefmt=\"none\")\n",
    "ax.set_xlim(0.0, 1000)\n",
    "ax.axvline(38, ls=\"--\", c=\"k\")\n",
    "ax.set_xlabel(r\"$f$ in $Hz$\")\n",
    "ax.set_ylabel(\"amplitude\")\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "**Congratulations! This completes the seventh and eighth exercise session.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nasty-giant",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}