notebooks/oat15_local_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fddbc3f0-63d1-4038-863e-6f6ad8302868",
   "metadata": {},
   "source": [
    "# Influence of weighting on the local prediction error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bcb9707-2bb1-4fa3-a9f6-50cb48b7137b",
   "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 as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "from flowtorch.analysis import DMD\n",
    "from utils import add_oat_patch\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/local_error\"\n",
    "\n",
    "makedirs(path, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a86a6d5-b45b-4c60-99f3-2e5d7b9753a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "x = vertices[\"x_small\"]\n",
    "z = vertices[\"z_small\"]\n",
    "del vertices\n",
    "start_at, end_at = 101, 501 # encloses 2 cycles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb17e6e2-1eb7-43a5-8407-4ef24d0e84f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "vel_x = pt.load(join(data, \"vel_x_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",
    "dm = pt.cat((vel_x, vel_z), dim=0)\n",
    "del vel_x, vel_z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "033982ea-ab1d-4599-aa80-e191f4912177",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = times[1] - times[0]\n",
    "rank = 65\n",
    "dmd = DMD(dm, dt, rank=rank, optimal=True)\n",
    "weights = area.sqrt().unsqueeze(-1).repeat((2, 1))\n",
    "dmdw = DMD(dm*weights, dt, rank=rank, optimal=True)\n",
    "print(dmd, dmdw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48f4c671-f9df-4924-bbfb-4754cf38cf36",
   "metadata": {},
   "outputs": [],
   "source": [
    "YH = (dmd.modes @ pt.diag(dmd.eigvals)) @ (pt.linalg.pinv(dmd.modes) @ dmd._dm[:, :-1].type(dmd.modes.dtype))\n",
    "p_err = (YH.real - dm[:, 1:]).abs().sum(dim=1)\n",
    "YH = (dmdw.modes @ pt.diag(dmdw.eigvals)) @ (pt.linalg.pinv(dmdw.modes) @ dmdw._dm[:, :-1].type(dmdw.modes.dtype))\n",
    "p_err_w = (YH.real / weights - dm[:, 1:]).abs().sum(dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69ff7cf2-efe0-4e1d-8e84-155e3a3be639",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_err.min(), p_err.max(), p_err_w.min(), p_err_w.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eda1d83-ba64-4991-af08-c7b898fe0132",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fmt(x, pos):\n",
    "    a, b = '{:.0e}'.format(x).split('e')\n",
    "    b = int(b)\n",
    "    return r'${} \\times 10^{{{}}}$'.format(a, b)\n",
    "\n",
    "CHORD = 0.15\n",
    "U_INF = 238.59\n",
    "n_points = area.shape[0]\n",
    "fig, axarr = plt.subplots(3, 1, figsize=(6, 4), sharex=True, sharey=True)\n",
    "vmin, vmax = -250, 250\n",
    "levels = pt.linspace(vmin, vmax, 100)\n",
    "cont_0 = axarr[0].tricontourf(x/CHORD, z/CHORD, p_err[:n_points]-p_err_w[:n_points], cmap=\"seismic\", levels=levels, vmin=vmin, vmax=vmax, extend=\"both\")\n",
    "cbar_0 = plt.colorbar(cont_0, ax=axarr[0], label=r\"$E_{\\mathrm{cum}} - E_{\\mathrm{cum}}^w$\", extend=\"both\")\n",
    "cbar_0.ax.get_yaxis().set_ticks([vmin, 0, vmax])\n",
    "vmin, vmax = 5e-6, 1e-3\n",
    "levels = pt.linspace(vmin, vmax, 100)\n",
    "cont_1 = axarr[1].tricontourf(x/CHORD, z/CHORD, area.sqrt(), levels=levels, vmin=vmin, vmax=vmax, extend=\"max\")\n",
    "cbar_1 = plt.colorbar(cont_1, ax=axarr[1], label=r\"$\\sqrt{A}$\", extend=\"max\", format=ticker.FuncFormatter(fmt))\n",
    "cbar_1.ax.get_yaxis().set_ticks([5e-6, 5e-4, 1e-3])\n",
    "dm_std = dm.std(dim=1)[:n_points] / U_INF\n",
    "vmin, vmax = 0.0, 0.9*dm_std.max()\n",
    "levels = pt.linspace(vmin, vmax, 100)\n",
    "cont_2 = axarr[2].tricontourf(x/CHORD, z/CHORD, dm_std, levels=levels, vmin=vmin, vmax=vmax, extend=\"max\")\n",
    "cbar_2 = plt.colorbar(cont_2, ax=axarr[2], label=r\"$\\sigma_{u_x}/U_\\infty$\", extend=\"max\", format=r\"$%1.2f$\")\n",
    "cbar_2.ax.get_yaxis().set_ticks([vmin, 0.5*(vmin+vmax), vmax])\n",
    "for ax in axarr:\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",
    "fig.add_subplot(111, frame_on=False)\n",
    "plt.tick_params(labelcolor=\"none\", bottom=False, left=False)\n",
    "plt.ylabel(r\"$\\tilde{z}$\", labelpad=-10)\n",
    "axarr[-1].set_xlabel(r\"$\\tilde{x}$\")\n",
    "axarr[0].text(2.3, 0.4, \"(a)\")\n",
    "axarr[1].text(2.3, 0.4, \"(b)\")\n",
    "axarr[2].text(2.3, 0.4, \"(c)\", c=\"w\")\n",
    "plt.savefig(join(path, \"projection_error_comparison.png\"), bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dda405b1-a9ff-4648-a09f-63c6ca69e49d",
   "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
}