notebooks/oat15_dmd_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9c8d467f-7284-4592-bff9-4d12b66ab030",
   "metadata": {},
   "source": [
    "# OAT15 DMD analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54322590-6865-4859-951a-a8d84007a83e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "from os import makedirs\n",
    "import torch as pt\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import FancyArrowPatch\n",
    "from flowtorch.analysis import DMD\n",
    "from utils import normalize_frequency, add_oat_patch, spatio_temporal_correlation\n",
    "\n",
    "plt.rcParams['figure.dpi'] = 160\n",
    "plt.rc('text', usetex=True)\n",
    "\n",
    "data = \"/home/andre/Development/naca0012_shock_buffet/run/oat15\"\n",
    "path = \"./output/oat15_analysis/dmd_analysis\"\n",
    "\n",
    "makedirs(path, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce7438e5-8822-48be-8efc-c11c7cfa276f",
   "metadata": {},
   "outputs": [],
   "source": [
    "CHORD = 0.15\n",
    "U_INF = 238.59\n",
    "\n",
    "times = pt.load(join(data, \"oat15_tandem_times.pt\"))[::20]\n",
    "vertices = pt.load(join(data, \"vertices_and_masks.pt\"))\n",
    "area = vertices[\"area_small\"]\n",
    "weight = area.sqrt().unsqueeze(-1)\n",
    "x = vertices[\"x_small\"] / CHORD\n",
    "z = vertices[\"z_small\"] / CHORD\n",
    "del vertices\n",
    "start_at, end_at = 101, 501 # encloses 2 cycles\n",
    "dt = times[1] - times[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30d0b70d-680d-44e6-b30a-041a4642d0f6",
   "metadata": {},
   "source": [
    "## Density - weighted and unweighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6093957-649f-4bf7-bfb3-0e9ae131af7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = pt.load(join(data, \"rho_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "dmd = DMD(dm, dt, optimal=True)\n",
    "dmdw = DMD(dm*weight, dt, optimal=True)\n",
    "print(dmd, dmdw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c775b2e-f512-413f-b3c2-05dc981f3ee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 2))\n",
    "\n",
    "freq = normalize_frequency(dmd.frequency, CHORD, U_INF)\n",
    "imp = dmd.integral_contribution[freq > 0.001]\n",
    "imp /= imp.max()\n",
    "freq = freq[freq > 0.001]\n",
    "freqw = normalize_frequency(dmdw.frequency, CHORD, U_INF)\n",
    "impw = dmdw.integral_contribution[freqw > 0.001]\n",
    "impw /= impw.max()\n",
    "freqw = freqw[freqw > 0.001]\n",
    "\n",
    "markerlinew, stemlinesw, baselinew = ax.stem(freqw.numpy(), impw.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C1o\", label=\"weighted\")\n",
    "plt.setp(stemlinesw, 'linewidth', 1.0)\n",
    "plt.setp(markerlinew, markersize=3)\n",
    "markerline, stemlines, baseline = ax.stem(freq.numpy(), imp.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C0x\", label=\"unweighted\")\n",
    "plt.setp(stemlines, 'linewidth', 1.0)\n",
    "plt.setp(markerline, markersize=3)\n",
    "ax.set_xlim(0, 20)\n",
    "ax.set_ylim(0, 1.05)\n",
    "ax.legend()\n",
    "ax.set_xlabel(r\"$\\tilde{f}$\")\n",
    "ax.set_ylabel(r\"$I/I_{max}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ae6cfcc-fa88-4784-8d87-5756ae8b4645",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_mode_symmetric(ax, mode, n_levels=60):\n",
    "    vmin, vmax = mode.min(), mode.max()\n",
    "    vmax = max(vmin.abs(), vmax.abs())\n",
    "    vmin = -vmax\n",
    "    levels = pt.linspace(vmin, vmax, n_levels)\n",
    "    ax.tricontourf(x, z, mode, levels=levels, vmin=vmin, vmax=vmax, cmap=\"seismic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cbc60b2-91b3-469f-aa89-578706499f2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_modes = 10\n",
    "top = dmd.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "topw = dmdw.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "fig, axarr = plt.subplots(n_modes, 2, figsize=(6, 1*n_modes), sharex=True, sharey=True)\n",
    "\n",
    "for i, (mwi, mi) in enumerate(zip(topw, top)):\n",
    "    add_mode_symmetric(axarr[i, 0], dmdw.modes[:, mwi].real/weight.squeeze())\n",
    "    add_mode_symmetric(axarr[i, 1], dmd.modes[:, mi].real)\n",
    "    axarr[i, 0].set_title(r\"$\\Re (\\rho_w)$\" + f\", {normalize_frequency(dmdw.frequency[mwi], CHORD, U_INF):2.2f}\")\n",
    "    axarr[i, 1].set_title(r\"$\\Re (\\rho)$\" + f\", {normalize_frequency(dmd.frequency[mi], CHORD, U_INF):2.2f}\")\n",
    "    for ax in axarr.flatten():\n",
    "        ax.set_aspect(\"equal\")\n",
    "        add_oat_patch(ax)\n",
    "        ax.set_ylim(-0.2, 0.5)\n",
    "        ax.set_xlim(-0.1, 2.5)\n",
    "\n",
    "for ax in axarr[:, 0]:\n",
    "    ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "for ax in axarr[-1, :]:\n",
    "    ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1fb2af8-f5c1-499b-943a-b593508212cf",
   "metadata": {},
   "source": [
    "## Pressure weighted and unweighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d1973f-d15e-4db8-87ca-f2dd503c7b15",
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = pt.load(join(data, \"p_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "dmd = DMD(dm, dt, optimal=True)\n",
    "dmdw = DMD(dm*weight, dt, optimal=True)\n",
    "print(dmd, dmdw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ae8fd3e-d58c-4857-9562-14c871640fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 2))\n",
    "\n",
    "freq = normalize_frequency(dmd.frequency, CHORD, U_INF)\n",
    "imp = dmd.integral_contribution[freq > 0.001]\n",
    "imp /= imp.max()\n",
    "freq = freq[freq > 0.001]\n",
    "freqw = normalize_frequency(dmdw.frequency, CHORD, U_INF)\n",
    "impw = dmdw.integral_contribution[freqw > 0.001]\n",
    "impw /= impw.max()\n",
    "freqw = freqw[freqw > 0.001]\n",
    "\n",
    "markerlinew, stemlinesw, baselinew = ax.stem(freqw.numpy(), impw.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C1o\", label=\"weighted\")\n",
    "plt.setp(stemlinesw, 'linewidth', 1.0)\n",
    "plt.setp(markerlinew, markersize=3)\n",
    "markerline, stemlines, baseline = ax.stem(freq.numpy(), imp.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C0x\", label=\"unweighted\")\n",
    "plt.setp(stemlines, 'linewidth', 1.0)\n",
    "plt.setp(markerline, markersize=3)\n",
    "ax.set_xlim(0, 20)\n",
    "ax.set_ylim(0, 1.05)\n",
    "ax.legend()\n",
    "ax.set_xlabel(r\"$\\tilde{f}$\")\n",
    "ax.set_ylabel(r\"$I/I_{max}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "066cbc5d-b68d-49b0-a954-5fffdec5c2f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_modes = 10\n",
    "top = dmd.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "topw = dmdw.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "fig, axarr = plt.subplots(n_modes, 2, figsize=(6, 1*n_modes), sharex=True, sharey=True)\n",
    "\n",
    "for i, (mwi, mi) in enumerate(zip(topw, top)):\n",
    "    add_mode_symmetric(axarr[i, 0], dmdw.modes[:, mwi].real/weight.squeeze())\n",
    "    add_mode_symmetric(axarr[i, 1], dmd.modes[:, mi].real)\n",
    "    axarr[i, 0].set_title(r\"$\\Re (p_w)$\" + f\", {normalize_frequency(dmdw.frequency[mwi], CHORD, U_INF):2.2f}\")\n",
    "    axarr[i, 1].set_title(r\"$\\Re (p)$\" + f\", {normalize_frequency(dmd.frequency[mi], CHORD, U_INF):2.2f}\")\n",
    "    for ax in axarr.flatten():\n",
    "        ax.set_aspect(\"equal\")\n",
    "        add_oat_patch(ax)\n",
    "        ax.set_ylim(-0.2, 0.5)\n",
    "        ax.set_xlim(-0.1, 2.5)\n",
    "\n",
    "for ax in axarr[:, 0]:\n",
    "    ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "for ax in axarr[-1, :]:\n",
    "    ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a804cd54-187c-41b2-bfbb-eeccfeec409f",
   "metadata": {},
   "source": [
    "## Velocity and local speed of sound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1edc9575-dc58-45b5-b201-8679821f1ada",
   "metadata": {},
   "outputs": [],
   "source": [
    "vel_x = pt.load(join(data, \"vel_x_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "vel_y = pt.load(join(data, \"vel_y_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "vel_z = pt.load(join(data, \"vel_z_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "ma = pt.load(join(data, \"ma_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "speed = (vel_x**2 + vel_y**2 + vel_z**2).sqrt()\n",
    "a_loc = speed / ma\n",
    "kappa = pt.tensor(1.4)\n",
    "scale = pt.sqrt(2.0 / (kappa * (kappa - 1.0)))\n",
    "dm_axz = pt.cat((vel_x, vel_z, a_loc*scale), dim=0)\n",
    "dm_xz = pt.cat((vel_x, vel_z))\n",
    "del vel_x, vel_y, vel_z, ma, a_loc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2fab33d-22a7-4213-a9ec-c368caa986ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "dmd_xz = DMD(dm_xz*weight.repeat((2, 1)), dt, optimal=True)\n",
    "dmd_axz = DMD(dm_axz*weight.repeat((3, 1)), dt, optimal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f5defb8-57a4-4b63-9a8d-1e3c2efe8dc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 2))\n",
    "\n",
    "freq_xz = normalize_frequency(dmd_xz.frequency, CHORD, U_INF)\n",
    "imp_xz = dmd_xz.integral_contribution[freq_xz > 0.001]\n",
    "imp_xz /= imp_xz.max()\n",
    "freq_xz = freq_xz[freq_xz > 0.001]\n",
    "freq_axz = normalize_frequency(dmd_axz.frequency, CHORD, U_INF)\n",
    "imp_axz = dmd_axz.integral_contribution[freq_axz > 0.001]\n",
    "imp_axz /= imp_axz.max()\n",
    "freq_axz = freq_axz[freq_axz > 0.001]\n",
    "\n",
    "markerline, stemlines, baseline = ax.stem(freq_xz.numpy(), imp_xz.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C1o\", label=r\"$\\mathbf{x}=(u_x, u_z)^T$\")\n",
    "plt.setp(stemlines, 'linewidth', 1.0)\n",
    "plt.setp(markerline, markersize=3)\n",
    "markerline, stemlines, baseline = ax.stem(freq_axz.numpy(), imp_axz.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C0x\", label=r\"$\\mathbf{x}=(u_x, u_z, \\sqrt{2/(\\gamma (\\gamma -1))} a)^T$\")\n",
    "plt.setp(stemlines, 'linewidth', 1.0)\n",
    "plt.setp(markerline, markersize=3)\n",
    "ax.set_xlim(0, 20)\n",
    "ax.set_ylim(0, 1.05)\n",
    "ax.legend()\n",
    "ax.set_xlabel(r\"$\\tilde{f}$\")\n",
    "ax.set_ylabel(r\"$I/I_{max}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18639b08-d207-460d-8c2d-bf14c8d86252",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_modes = 10\n",
    "n_points = weight.shape[0]\n",
    "top_xz = dmd_xz.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "top_axz = dmd_axz.top_modes(n_modes, integral=True, f_min=0.001)\n",
    "fig, axarr = plt.subplots(n_modes, 2, figsize=(6, 1*n_modes), sharex=True, sharey=True)\n",
    "\n",
    "for i, (mi_xz, mi_axz) in enumerate(zip(top_xz, top_axz)):\n",
    "    add_mode_symmetric(axarr[i, 0], dmd_xz.modes[:n_points, mi_xz].real/weight.squeeze())\n",
    "    add_mode_symmetric(axarr[i, 1], dmd_axz.modes[:n_points, mi_axz].real/weight.squeeze())\n",
    "    axarr[i, 0].set_title(r\"$\\Re (u_x)$\" + f\", {normalize_frequency(dmd_xz.frequency[mi_xz], CHORD, U_INF):2.2f}\")\n",
    "    axarr[i, 1].set_title(r\"$\\Re (u_x^a)$\" + f\", {normalize_frequency(dmd_axz.frequency[mi_axz], CHORD, U_INF):2.2f}\")\n",
    "    for ax in axarr.flatten():\n",
    "        ax.set_aspect(\"equal\")\n",
    "        add_oat_patch(ax)\n",
    "        ax.set_ylim(-0.2, 0.5)\n",
    "        ax.set_xlim(-0.1, 2.5)\n",
    "\n",
    "for ax in axarr[:, 0]:\n",
    "    ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "for ax in axarr[-1, :]:\n",
    "    ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f6f3a8-0bb4-4795-9690-452974d38bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 2))\n",
    "\n",
    "freq_xz = normalize_frequency(dmd_xz.frequency, CHORD, U_INF)\n",
    "imp_xz = dmd_xz.integral_contribution[freq_xz > 0.001]\n",
    "imp_xz /= imp_xz.max()\n",
    "freq_xz = freq_xz[freq_xz > 0.001]\n",
    "\n",
    "markerline, stemlines, baseline = ax.stem(freq_xz.numpy(), imp_xz.numpy(), basefmt=\"none\", linefmt=f\"k-\", markerfmt=f\"C1o\", label=\"opt. DMD, int.\")\n",
    "plt.setp(stemlines, 'linewidth', 1.0)\n",
    "plt.setp(markerline, markersize=4)\n",
    "fb = normalize_frequency(118.5, CHORD, U_INF)\n",
    "ax.axvline(fb, ls=\":\", c=\"k\", label=r\"$f_b$ (lift)\")\n",
    "ax.set_title(r\"$\\mathbf{x}=(u_x, u_z)^T$\")\n",
    "ax.set_xlim(0, 20)\n",
    "ax.set_ylim(0, 1.05)\n",
    "ax.set_xlabel(r\"$\\tilde{f}$\")\n",
    "ax.set_ylabel(r\"$I/I_{\\mathrm{max}}$\")\n",
    "ax.legend()\n",
    "plt.savefig(join(path, \"freq_vel_xz.pdf\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f331c2e5-2387-4229-9958-1ef83f08b097",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_modes = 6\n",
    "n_points = weight.shape[0]\n",
    "top_xz = dmd_xz.top_modes(1, integral=True)\n",
    "top_xz = pt.cat((top_xz, dmd_xz.top_modes(3, integral=True, f_min=0.001, f_max=400)))\n",
    "top_xz = pt.cat((top_xz, dmd_xz.top_modes(2, integral=True, f_min=400)))\n",
    "fig, axarr = plt.subplots(n_modes, 2, figsize=(6, 1.1*n_modes), sharex=True, sharey=True)\n",
    "\n",
    "for i, mi_xz in enumerate(top_xz):\n",
    "    add_mode_symmetric(axarr[i, 0], dmd_xz.modes[:n_points, mi_xz].real/weight.squeeze(), 100)\n",
    "    add_mode_symmetric(axarr[i, 1], dmd_xz.modes[n_points:, mi_xz].real/weight.squeeze(), 100)\n",
    "    axarr[i, 0].set_title(r\"$\\Re (u_x)$\" + r\", $\\tilde{f} = \" + \"{:2.2f}$\".format(normalize_frequency(dmd_xz.frequency[mi_xz], CHORD, U_INF)))\n",
    "    axarr[i, 1].set_title(r\"$\\Re (u_z)$\" + r\", $\\tilde{f} = \" + \"{:2.2f}$\".format(normalize_frequency(dmd_xz.frequency[mi_xz], CHORD, U_INF)))\n",
    "    for ax in axarr.flatten():\n",
    "        ax.set_aspect(\"equal\")\n",
    "        add_oat_patch(ax)\n",
    "        ax.set_ylim(-0.2, 0.5)\n",
    "        ax.set_xlim(-0.1, 2.5)\n",
    "\n",
    "for ax in axarr[-1, :]:\n",
    "    ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "    \n",
    "fig.add_subplot(111, frame_on=False)\n",
    "plt.tick_params(labelcolor=\"none\", bottom=False, left=False)\n",
    "plt.ylabel(r\"$\\tilde{z}$\")\n",
    "plt.savefig(join(path, \"vel_xz_modes.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05dea1e7-bd69-4d29-a8fa-c2f4cc18012d",
   "metadata": {},
   "source": [
    "## Propagation speeds - vortex shedding mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "101ecdf9-d1e2-44c4-b463-69f86e52d573",
   "metadata": {},
   "outputs": [],
   "source": [
    "top_xz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70311c1c-1d9e-4f1d-9295-1fb3a22e48d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "top_xz_conj = [top_xz[0]]\n",
    "for mi in top_xz[1:]:\n",
    "    top_xz_conj.extend((dmd_xz.frequency.abs() - dmd_xz.frequency[mi].abs()).abs().topk(2, largest=False).indices)\n",
    "top_xz_conj = pt.tensor(top_xz_conj, dtype=pt.int64)\n",
    "top_xz_conj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b891129a-e5f2-4fce-9912-9f4f95692e23",
   "metadata": {},
   "outputs": [],
   "source": [
    "top_rec_err = (dmd_xz.partial_reconstruction(top_xz_conj) - dmd_xz._dm).norm()\n",
    "max_rec_err = (dmd_xz.partial_reconstruction([top_xz[0]]) - dmd_xz._dm).norm()\n",
    "min_rec_err = dmd_xz.reconstruction_error.norm()\n",
    "top_rec_err, max_rec_err, min_rec_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa250e9-67bd-4b3e-9fd6-f1964ac73a95",
   "metadata": {},
   "outputs": [],
   "source": [
    "(top_rec_err - min_rec_err) / (max_rec_err - min_rec_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d04af13e-a478-435c-b434-2108697db0fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n",
    "mi = top_xz[4]\n",
    "phase_ux = dmd_xz.modes[:n_points, mi].angle()\n",
    "vmin, vmax= -180, 180\n",
    "levels = pt.linspace(vmin, vmax, 61)\n",
    "tri_ux = ax.tricontourf(x, z, phase_ux*180/np.pi, levels=levels, cmap=\"seismic\", vmin=vmin, vmax=vmax)\n",
    "bar = plt.colorbar(tri_ux, ax=ax, label=r\"$\\mathrm{tan}^{-1}\\left(\\Im (u_x) / \\Re (u_x)\\right)$\", shrink=0.8)\n",
    "ax.set_aspect(\"equal\")\n",
    "add_oat_patch(ax)\n",
    "ax.set_ylim(-0.5, 0.75)\n",
    "ax.set_xlim(-0.1, 2.5)\n",
    "ax.add_patch(\n",
    "    FancyArrowPatch((1.0, 0.5), (0.5, 0.5), mutation_scale=20, fc=\"w\")\n",
    ")\n",
    "ax.add_patch(\n",
    "    FancyArrowPatch((1.0, -0.3), (0.5, -0.3), mutation_scale=20, fc=\"w\")\n",
    ")\n",
    "ax.add_patch(\n",
    "    FancyArrowPatch((1.5, 0.0), (2.0, 0.0), mutation_scale=20, fc=\"w\")\n",
    ")\n",
    "ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "ax.set_title(r\"phase shift, $\\tilde{{f}}={:2.2f}$\".format(normalize_frequency(dmd_xz.frequency[mi], CHORD, U_INF)))\n",
    "\n",
    "plt.savefig(join(path, \"phase_shift.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ce5504-285c-4376-8472-8574e0c52cdc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 1.5))\n",
    "\n",
    "phase_ux = dmd_xz.modes[:n_points, mi].angle()\n",
    "vmin, vmax= -180, 180\n",
    "levels = pt.linspace(vmin, vmax, 61)\n",
    "tri_ux = ax.tricontourf(x, z, phase_ux*180/np.pi, levels=levels, cmap=\"seismic\", vmin=vmin, vmax=vmax)\n",
    "bar = plt.colorbar(tri_ux, ax=ax, label=r\"$\\mathrm{tan}^{-1}\\left(\\Im (u_x) / \\Re (u_x)\\right)$\", shrink=0.8)\n",
    "ax.set_aspect(\"equal\")\n",
    "add_oat_patch(ax)\n",
    "ax.set_ylim(0.0, 0.15)\n",
    "ax.set_xlim(0.5, 1.1)\n",
    "ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "ax.set_title(r\"phase shift, $\\tilde{{f}}={:2.2f}$\".format(normalize_frequency(dmd_xz.frequency[mi], CHORD, U_INF)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4fd8ee2-cb99-4b28-8961-34b55917cf45",
   "metadata": {},
   "outputs": [],
   "source": [
    "np_line = 45\n",
    "line = pt.stack((\n",
    "    pt.linspace(0.525, 0.72, np_line),\n",
    "    pt.linspace(0.3, 0.3, np_line),\n",
    ")).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0792f6d-2da0-4259-9f5a-e280af65631e",
   "metadata": {},
   "outputs": [],
   "source": [
    "phase_l = []\n",
    "for pi in range(line.shape[0]):\n",
    "    closest = (pt.stack((x, z), dim=1) - line[pi]).norm(dim=1).argmin()\n",
    "    phase_l.append(phase_ux[closest])\n",
    "phase_l = pt.stack(phase_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cfc1f10-006d-4ba3-81fa-9c0c81894054",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 2))\n",
    "ax.scatter(line[:, 0], phase_l, marker=\"x\", c=\"k\", s=10, label=r\"line sample, $\\tilde{z}=0.3$\")\n",
    "low_1, up_1 = 0.525, 0.61\n",
    "l1_range = pt.logical_and(line[:, 0] > low_1, line[:, 0] < up_1)\n",
    "dm = pt.cat((phase_l[l1_range].unsqueeze(-1), line[l1_range, 0].unsqueeze(-1)), dim=1).T\n",
    "dm_mean = dm.mean(dim=1)\n",
    "U, s, VH = pt.linalg.svd(dm - dm_mean.unsqueeze(-1), full_matrices=False)\n",
    "dphi_dx = U[0, 0] / U[1, 0]\n",
    "u_rel_1 = -2*np.pi*dmd_xz.frequency[mi] / dphi_dx / U_INF * CHORD\n",
    "phi_1_rec = dm_mean[0] + dphi_dx * (line[l1_range, 0] - dm_mean[1])\n",
    "ax.plot(line[l1_range, 0], phi_1_rec, ls=\"--\", c=\"C0\", label=\"linear fit\")\n",
    "ax.text(0.55, -0.5, r\"$U_p=\" + f\"{u_rel_1:2.2f}\" + \"U_\\infty$\", rotation=35)\n",
    "\n",
    "low_2, up_2 = 0.613, 0.703\n",
    "l2_range = pt.logical_and(line[:, 0] > low_2, line[:, 0] < up_2)\n",
    "dm = pt.cat((phase_l[l2_range].unsqueeze(-1), line[l2_range, 0].unsqueeze(-1)), dim=1).T\n",
    "dm_mean = dm.mean(dim=1)\n",
    "U, s, VH = pt.linalg.svd(dm - dm_mean.unsqueeze(-1), full_matrices=False)\n",
    "dphi_dx = U[0, 0] / U[1, 0]\n",
    "u_rel_2 = -2*np.pi*dmd_xz.frequency[mi] / dphi_dx / U_INF * CHORD\n",
    "phi_2_rec = dm_mean[0] + dphi_dx * (line[l2_range, 0] - dm_mean[1])\n",
    "ax.plot(line[l2_range, 0], phi_2_rec, ls=\"--\", c=\"C0\")\n",
    "ax.text(0.635, -0.5, r\"$U_p=\" + f\"{u_rel_2:2.2f}\" + \"U_\\infty$\", rotation=35)\n",
    "ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "ax.set_ylabel(r\"$\\mathrm{tan}^{-1}\\left(\\Im (u_x) / \\Re (u_x)\\right)$\")\n",
    "ax.legend()\n",
    "#ax.axvline(0.5275, ls=\":\", c=\"k\")\n",
    "#ax.axvline(0.6115, ls=\":\", c=\"k\")\n",
    "#ax.text(0.55, -2, r\"$\\Delta x \\approx {:2.3f}c$\".format(0.6115-0.5275))\n",
    "plt.savefig(join(path, \"phase_speed.pdf\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ecfa1a4-bb4e-4b41-a42f-2935a8e30966",
   "metadata": {},
   "source": [
    "## Propagation speeds - buffet mode (based on pressure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2450646-842c-4431-bac4-9fe64efdab18",
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = pt.load(join(data, \"p_small_every10.pt\"))[:, start_at:end_at:2]\n",
    "dmdw = DMD(dm*weight, dt, optimal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29dba7ce-52f7-476e-9192-440c8e2bff37",
   "metadata": {},
   "outputs": [],
   "source": [
    "topw = dmdw.top_modes(10, integral=True, f_min=0.001)\n",
    "topw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3980f292-3bfb-494f-b373-078a27da5155",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalize_frequency(dmdw.frequency[topw], CHORD, U_INF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "854eb276-f62d-4218-806d-9969283f60eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n",
    "mi = topw[0]\n",
    "phase_p = dmdw.modes[:, mi].angle()\n",
    "vmin, vmax= -180, 180\n",
    "levels = pt.linspace(vmin, vmax, 61)\n",
    "tri_p = ax.tricontourf(x, z, phase_p*180/np.pi, levels=levels, cmap=\"seismic\", vmin=vmin, vmax=vmax)\n",
    "bar = plt.colorbar(tri_p, ax=ax, label=r\"$\\mathrm{tan}^{-1}\\left(\\Im (p) / \\Re (p)\\right)$\", shrink=0.8)\n",
    "ax.set_aspect(\"equal\")\n",
    "add_oat_patch(ax)\n",
    "ax.set_ylim(-0.5, 0.75)\n",
    "ax.set_xlim(-0.1, 2.5)\n",
    "ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "ax.set_title(r\"phase shift, $\\tilde{{f}}={:2.2f}$\".format(normalize_frequency(dmdw.frequency[mi], CHORD, U_INF)))\n",
    "\n",
    "plt.savefig(join(path, \"phase_shift_buffet.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6b3ed73-330a-410b-bdf0-2778cd9e960b",
   "metadata": {},
   "outputs": [],
   "source": [
    "np_line_1 = 20\n",
    "line_1 = pt.stack((\n",
    "    pt.linspace(0.6, 0.9, np_line_1),\n",
    "    pt.linspace(0.078, 0.033, np_line_1),\n",
    ")).T\n",
    "np_line_2 = 10\n",
    "line_2 = pt.stack((\n",
    "    pt.linspace(0.88, 1.0, np_line_2),\n",
    "    pt.linspace(-0.008, -0.015, np_line_2),\n",
    ")).T\n",
    "np_line_3 = 15\n",
    "line_3 = pt.stack((\n",
    "    pt.linspace(0.3, 0.8, np_line_3),\n",
    "    pt.linspace(-0.08, -0.015, np_line_3),\n",
    ")).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eef30c12-180f-4b46-930d-7b1d1f4a83fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n",
    "\n",
    "vmin, vmax= -180, 180\n",
    "levels = pt.linspace(vmin, vmax, 61)\n",
    "tri_p = ax.tricontourf(x, z, phase_p*180/np.pi, levels=levels, cmap=\"seismic\", vmin=vmin, vmax=vmax)\n",
    "bar = plt.colorbar(tri_p, ax=ax, label=r\"$\\mathrm{tan}^{-1}\\left(\\Im (p) / \\Re (p)\\right)$\", shrink=0.8)\n",
    "ax.set_aspect(\"equal\")\n",
    "add_oat_patch(ax)\n",
    "ax.set_ylim(-0.15, 0.4)\n",
    "ax.set_xlim(0.0, 1.1)\n",
    "ax.plot(line_1[:, 0], line_1[:, 1], lw=1, c=\"k\", ls=\"--\")\n",
    "ax.plot(line_2[:, 0], line_2[:, 1], lw=1, c=\"k\", ls=\"--\")\n",
    "ax.plot(line_3[:, 0], line_3[:, 1], lw=1, c=\"k\", ls=\"--\")\n",
    "ax.set_xlabel(r\"$\\tilde{x}$\")\n",
    "ax.set_ylabel(r\"$\\tilde{z}$\")\n",
    "ax.text(0.45, -0.1, r\"$-0.25 U_\\infty$\", rotation=7)\n",
    "ax.text(0.7, 0.055, r\"$0.070 U_\\infty$\", rotation=-10)\n",
    "ax.text(0.85, -0.06, r\"$-0.17 U_\\infty$\", rotation=-5)\n",
    "ax.set_title(r\"phase shift, $\\tilde{{f}}={:2.2f}$\".format(normalize_frequency(dmdw.frequency[mi], CHORD, U_INF)))\n",
    "plt.savefig(join(path, \"phase_shift_buffet_enlarged.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc96d507-10c4-490c-bbe9-0f27c914bb1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "phase_l1 = []\n",
    "for pi in range(line_1.shape[0]):\n",
    "    closest = (pt.stack((x, z), dim=1) - line_1[pi]).norm(dim=1).argmin()\n",
    "    phase_l1.append(phase_p[closest])\n",
    "phase_l1 = pt.stack(phase_l1)\n",
    "phase_l2 = []\n",
    "for pi in range(line_2.shape[0]):\n",
    "    closest = (pt.stack((x, z), dim=1) - line_2[pi]).norm(dim=1).argmin()\n",
    "    phase_l2.append(phase_p[closest])\n",
    "phase_l2 = pt.stack(phase_l2)\n",
    "phase_l3 = []\n",
    "for pi in range(line_3.shape[0]):\n",
    "    closest = (pt.stack((x, z), dim=1) - line_3[pi]).norm(dim=1).argmin()\n",
    "    phase_l3.append(phase_p[closest])\n",
    "phase_l3 = pt.stack(phase_l3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "393c737a-87d1-4601-b974-3f77f90ae8a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axarr = plt.subplots(1, 3, figsize=(6, 2), sharey=True, sharex=False)\n",
    "\n",
    "# line 1 - boundary layer suction side\n",
    "ds_1 = (line_1 - line_1[0, :]).norm(dim=1)\n",
    "axarr[0].scatter(ds_1, phase_l1, marker=\"x\", c=\"k\", s=10)\n",
    "low_1, up_1 = 0.0, 0.4\n",
    "l1_range = pt.logical_and(ds_1 > low_1, ds_1 < up_1)\n",
    "dm = pt.cat((phase_l1[l1_range].unsqueeze(-1), line_1[l1_range, 0].unsqueeze(-1)), dim=1).T\n",
    "dm_mean = dm.mean(dim=1)\n",
    "U, s, VH = pt.linalg.svd(dm - dm_mean.unsqueeze(-1), full_matrices=False)\n",
    "dphi_ds_1 = U[0, 0] / U[1, 0]\n",
    "u_rel_1 = -2*np.pi*dmdw.frequency[mi] / dphi_ds_1 / U_INF * CHORD\n",
    "phi_1_rec = dm_mean[0] + dphi_ds_1 * (line_1[l1_range, 0] - dm_mean[1])\n",
    "axarr[0].plot(ds_1[l1_range], phi_1_rec, ls=\"--\", c=\"C0\", label=\"linear fit\")\n",
    "axarr[0].set_title(r\"$U_p=\" + f\"{u_rel_1:2.3f}\" + \"U_\\infty$\")\n",
    "\n",
    "# line 2 - pressure side towards trailing edge\n",
    "ds_2 = (line_2 - line_2[0, :]).norm(dim=1)\n",
    "axarr[1].scatter(ds_2, phase_l2, marker=\"x\", c=\"k\", s=10)\n",
    "low_2, up_2 = 0.0, 0.15\n",
    "l2_range = pt.logical_and(ds_2 > low_2, ds_2 < up_2)\n",
    "dm = pt.cat((phase_l2[l2_range].unsqueeze(-1), line_2[l2_range, 0].unsqueeze(-1)), dim=1).T\n",
    "dm_mean = dm.mean(dim=1)\n",
    "U, s, VH = pt.linalg.svd(dm - dm_mean.unsqueeze(-1), full_matrices=False)\n",
    "dphi_ds_2 = U[0, 0] / U[1, 0]\n",
    "u_rel_2 = -2*np.pi*dmdw.frequency[mi] / dphi_ds_2 / U_INF * CHORD\n",
    "phi_2_rec = dm_mean[0] + dphi_ds_2 * (line_2[l2_range, 0] - dm_mean[1])\n",
    "axarr[1].plot(ds_2[l2_range], phi_2_rec, ls=\"--\", c=\"C0\", label=\"linear fit\")\n",
    "axarr[1].set_title(r\"$U_p=\" + f\"{u_rel_2:2.2f}\" + \"U_\\infty$\")\n",
    "\n",
    "# line 3 - pressure side center\n",
    "ds_3 = (line_3 - line_3[0, :]).norm(dim=1)\n",
    "axarr[2].scatter(ds_3, phase_l3, marker=\"x\", c=\"k\", s=10)\n",
    "low_3, up_3 = 0.0, 0.5\n",
    "l3_range = pt.logical_and(ds_3 > low_3, ds_3 < up_3)\n",
    "dm = pt.cat((phase_l3[l3_range].unsqueeze(-1), line_3[l3_range, 0].unsqueeze(-1)), dim=1).T\n",
    "dm_mean = dm.mean(dim=1)\n",
    "U, s, VH = pt.linalg.svd(dm - dm_mean.unsqueeze(-1), full_matrices=False)\n",
    "dphi_ds_3 = U[0, 0] / U[1, 0]\n",
    "u_rel_3 = -2*np.pi*dmdw.frequency[mi] / dphi_ds_3 / U_INF * CHORD\n",
    "phi_3_rec = dm_mean[0] + dphi_ds_3 * (line_3[l3_range, 0] - dm_mean[1])\n",
    "axarr[2].plot(ds_3[l3_range], phi_3_rec, ls=\"--\", c=\"C0\", label=\"linear fit\")\n",
    "axarr[2].set_title(r\"$U_p=\" + f\"{u_rel_3:2.2f}\" + \"U_\\infty$\")\n",
    "for ax in axarr.flatten():\n",
    "    ax.set_xlabel(r\"$\\tilde{s}$\")\n",
    "axarr[0].set_ylabel(r\"$\\mathrm{tan}^{-1}\\left(\\Im (p) / \\Re (p)\\right)$\")\n",
    "plt.savefig(join(path, \"phase_speed_buffet.pdf\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54589c72-4e15-447e-b44c-5eec0f54669c",
   "metadata": {},
   "source": [
    "## UDMD buffet mode - rank dependency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4a23b39-c9aa-43a9-9f9c-3b27f6dd4274",
   "metadata": {},
   "outputs": [],
   "source": [
    "ranks = [25, 50, 100]\n",
    "n_points = weight.shape[0]\n",
    "top_mode_dmd, top_mode_udmd = [], []\n",
    "for r in ranks:\n",
    "    dmd = DMD(dm_xz*weight.repeat((2, 1)), dt, rank=r, optimal=True)\n",
    "    top = dmd.top_modes(1, integral=True, f_min=1)\n",
    "    top_mode_dmd.append(dmd.modes[:n_points, top[0]].real)\n",
    "    dmd = DMD(dm_xz*weight.repeat((2, 1)), dt, rank=r, optimal=True, unitary=True)\n",
    "    top = dmd.top_modes(1, integral=True, f_min=1)\n",
    "    top_mode_udmd.append(dmd.modes[:n_points, top[0]].real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a32c1770-76b0-4741-800b-b4f5c3305bc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axarr = plt.subplots(2, 3, figsize=(6, 1.5), sharex=True, sharey=True)\n",
    "\n",
    "for i, r in enumerate(ranks):\n",
    "    add_mode_symmetric(axarr[0, i], top_mode_udmd[i] / weight.squeeze())\n",
    "    add_mode_symmetric(axarr[1, i], top_mode_dmd[i] / weight.squeeze())\n",
    "    axarr[0, i].set_title(r\"$\\Re (u_x)$, $r={:d}$\".format(r))\n",
    "    axarr[1, i].set_xlabel(r\"$\\tilde{x}$\")\n",
    "    for ax in axarr[:, i]:\n",
    "        ax.set_aspect(\"equal\")\n",
    "        add_oat_patch(ax)\n",
    "        ax.set_ylim(-0.2, 0.5)\n",
    "        ax.set_xlim(-0.1, 2.5)\n",
    "\n",
    "axarr[0, 0].set_ylabel(r\"$\\tilde{z}$\")\n",
    "axarr[1, 0].set_ylabel(r\"$\\tilde{z}$\")\n",
    "axarr[0, 2].text(2.6, 0.2, \"UDMD\", rotation=90, va=\"center\")\n",
    "axarr[1, 2].text(2.6, 0.2, \"DMD\", rotation=90, va=\"center\")\n",
    "plt.savefig(join(path, \"top_mode_udmd_dmd.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2385c190-a89c-40a3-9c75-9204a7c0c489",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}