analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Specifying the paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_path = Path(\n",
    "    r\"\\\\wsl.localhost\\Ubuntu-24.04\\home\\user\\OF\\attempt_2\\results\\csv\\results.csv\"\n",
    ")\n",
    "figures_path = Path(\n",
    "    r\"\\\\wsl.localhost\\Ubuntu-24.04\\home\\user\\OF\\attempt_2\\results/figures\"\n",
    ")\n",
    "figures_path.mkdir(exist_ok=True, parents=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "force_coeff_df = pd.read_csv(csv_path, parse_dates=[\"datetime\"])\n",
    "force_coeff_df[\"cl_cd\"] = force_coeff_df[\"cl\"] / force_coeff_df[\"cd\"]\n",
    "force_coeff_df[\"cl_cd_abs\"] = (force_coeff_df[\"cl\"] / force_coeff_df[\"cd\"]).abs()\n",
    "force_coeff_df[\"failed\"] = (\n",
    "    ~force_coeff_df[[\"no_clipping\", \"block_mesh\", \"check_mesh\", \"simple\"]]\n",
    ").any(axis=1)\n",
    "\n",
    "df_filtered = force_coeff_df\n",
    "\n",
    "df_filtered = force_coeff_df.dropna(subset=[\"cl\", \"cd\"]).copy()\n",
    "# df_filtered = df_filtered[\n",
    "#     df_filtered[\"run_name\"] == \"5_degree_AoA_fixed_nu_tilda_reduced_yplus\"\n",
    "# ]\n",
    "# df_filtered = df_filtered[\n",
    "#     df_filtered[\"cl_cd_abs\"] < 200\n",
    "# ]\n",
    "df_filtered[\"moving_avg\"] = df_filtered[\"cl_cd_abs\"].rolling(window=30).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filtered.sort_values(by=\"cl_cd_abs\", ascending=False)\n",
    "# len(df_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Progress plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.scatter(\n",
    "    df_filtered[\"datetime\"],\n",
    "    df_filtered[\"cl_cd_abs\"],\n",
    "    label=\"Lift-drag ratio (Cl/Cd)\",\n",
    "    marker=\"o\",\n",
    "    color=\"tab:blue\",\n",
    "    alpha=0.25,\n",
    ")\n",
    "plt.plot(\n",
    "    df_filtered[\"datetime\"],\n",
    "    df_filtered[\"cl_cd_abs\"].cummax(),\n",
    "    label=\"Cumulative max Cl/Cd\",\n",
    "    linestyle=\"--\",\n",
    "    color=\"tab:green\",\n",
    ")\n",
    "plt.scatter(\n",
    "    force_coeff_df[force_coeff_df[\"failed\"]][\"datetime\"],\n",
    "    -force_coeff_df[force_coeff_df[\"failed\"]][\"failed\"],\n",
    "    s=10,\n",
    "    color=\"tab:red\",\n",
    "    alpha=0.1,\n",
    "    label=\"Failed runs\",\n",
    ")\n",
    "\n",
    "plt.plot(\n",
    "    df_filtered[\"datetime\"],\n",
    "    df_filtered[\"moving_avg\"],\n",
    "    label=\"Moving average\",\n",
    "    linestyle=\"-.\",\n",
    "    color=\"tab:orange\",\n",
    ")\n",
    "\n",
    "plt.title(f\"Lift-drag ratio over time - {len(force_coeff_df)} simulations\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"$C_l/C_d$ ratio\")\n",
    "\n",
    "# Use a locator to reduce the density of x-ticks\n",
    "# locator = mdates.HourLocator(interval=1)\n",
    "# formatter = mdates.DateFormatter(\"%H:%M\")\n",
    "# plt.gca().xaxis.set_major_locator(locator)\n",
    "# plt.gca().xaxis.set_major_formatter(formatter)\n",
    "\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\n",
    "    \"results/figures/27122024 - OpenFOAM - Lift-drag ratio over time.png\", dpi=300\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Attempts at model reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_subset = df_filtered[\n",
    "    (df_filtered[\"cl_cd_abs\"] > 5.0) & (df_filtered[\"cl_cd_abs\"] < 1005.0)\n",
    "]\n",
    "\n",
    "X = df_subset[[\"x0\", \"x1\", \"x2\", \"x3\", \"x4\", \"x5\"]]\n",
    "y = df_subset[\"cl_cd\"]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.15)\n",
    "\n",
    "rf = RandomForestRegressor(n_estimators=100, random_state=42)\n",
    "\n",
    "# Train the model\n",
    "rf.fit(X_train, y_train)\n",
    "y_pred = rf.predict(X_test)\n",
    "\n",
    "mse = mean_squared_error(y_test, y_pred)\n",
    "mae = mean_absolute_error(y_test, y_pred)\n",
    "percentage_error = np.mean(abs((y_test - y_pred) / y_test)) * 100\n",
    "\n",
    "print(f\"Mean Squared Error: {mse}\")\n",
    "print(f\"Mean Absolute Error: {mae}\")\n",
    "print(f\"Mean Percentage Error: {percentage_error:.2f}%\")\n",
    "\n",
    "# Dummy method: Mean baseline\n",
    "y_dummy_pred = [y_train.mean()] * len(y_test)\n",
    "mse_dummy = mean_squared_error(y_test, y_dummy_pred)\n",
    "mae_dummy = mean_absolute_error(y_test, y_dummy_pred)\n",
    "percentage_error_dummy = np.mean(abs((y_test - y_dummy_pred) / y_test)) * 100\n",
    "\n",
    "print(\"\\nDummy Model Evaluation:\")\n",
    "print(f\"Mean Squared Error (Dummy): {mse_dummy}\")\n",
    "print(f\"Mean Absolute Error (Dummy): {mae_dummy}\")\n",
    "print(f\"Mean Percentage Error (Dummy): {percentage_error_dummy:.2f}%\")\n",
    "\n",
    "# plt.hist((y_test - y_pred) / y_test * 100, bins=100)\n",
    "plt.hist(abs((y_test - y_pred)), bins=30)\n",
    "plt.xlabel(\"Error (MAE)\")\n",
    "plt.ylabel(\"Bin count\")\n",
    "plt.title(\"Model reduction with random forest\")\n",
    "plt.grid()\n",
    "\n",
    "plt.savefig(\n",
    "    \"results/figures/27122024 - OpenFOAM - random forest model reduction.png\", dpi=300\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Predicted vs. actual values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option(\"display.max_rows\", None)\n",
    "model_reduction_df = pd.DataFrame({\"pred\": y_pred, \"test\": y_test.values}).sort_values(\n",
    "    by=\"test\", ascending=False\n",
    ")\n",
    "model_reduction_df[\"error\"] = model_reduction_df[\"test\"] - model_reduction_df[\"pred\"]\n",
    "model_reduction_df[\"percent_error\"] = (\n",
    "    (model_reduction_df[\"test\"] - model_reduction_df[\"pred\"])\n",
    "    / model_reduction_df[\"test\"]\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Forces over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "force_coeff_df = pd.read_csv(\n",
    "    r\"\\\\wsl.localhost\\Ubuntu-24.04\\home\\user\\OF\\attempt_2\\custom_runs/aoa\\41.05_degree_AoA\\postProcessing\\forceCoeffs\\0\\coefficient.dat\",\n",
    "    skiprows=12,\n",
    "    sep=\"\\t\",\n",
    ")\n",
    "force_coeff_df.columns = [column_name.strip() for column_name in force_coeff_df.columns]\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "# Plot Cl and Cd on the first y-axis\n",
    "(line1,) = ax1.plot(\n",
    "    force_coeff_df[\"# Time\"], force_coeff_df[\"Cl\"], label=\"$C_l$\", color=\"tab:blue\"\n",
    ")\n",
    "(line2,) = ax1.plot(\n",
    "    force_coeff_df[\"# Time\"], force_coeff_df[\"Cd\"], label=\"$C_d$\", color=\"tab:green\"\n",
    ")\n",
    "ax1.set_xlabel(\"Time\")\n",
    "ax1.set_ylabel(\"$C_l$, $C_d$\")\n",
    "# ax1.set_ylim([-0.0005, 0.0005])\n",
    "\n",
    "# Create a second y-axis for Cl/Cd\n",
    "ax2 = ax1.twinx()\n",
    "(line3,) = ax2.plot(\n",
    "    force_coeff_df[\"# Time\"],\n",
    "    force_coeff_df[\"Cl\"] / force_coeff_df[\"Cd\"],\n",
    "    label=\"$C_l/C_d$\",\n",
    "    color=\"tab:red\",\n",
    ")\n",
    "ax2.set_ylabel(\"$C_l/C_d$\")\n",
    "ax2.set_ylim([0, 50])\n",
    "\n",
    "# Combine legends\n",
    "lines = [line1, line2, line3]\n",
    "labels = [line.get_label() for line in lines]\n",
    "ax1.legend(lines, labels, loc=\"upper right\")\n",
    "\n",
    "plt.title(\"Forces over time\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\n",
    "    \"results/figures/01012025 - OpenFOAM - oscillatory solution to lift-drag ratio - 41.05 deg AoA.png\",\n",
    "    dpi=300,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lift-drag vs angle-of-attack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aoa_df = pd.read_csv(\"results/csv/aoa_results.csv\")\n",
    "\n",
    "aoa_df[\"aoa\"] = aoa_df[\"run_name\"].str.extract(r\"(\\d+\\.\\d+)_degree_AoA\").astype(float)\n",
    "aoa_df = aoa_df[aoa_df[\"cl\"] >= 0]\n",
    "\n",
    "aoa_df = aoa_df.sort_values(by=\"aoa\")\n",
    "\n",
    "plt.plot(aoa_df[\"aoa\"], aoa_df[\"cl\"] / aoa_df[\"cd\"], \"-x\")\n",
    "plt.xlabel(\"Angle-of-attack\")\n",
    "plt.ylabel(\"$C_l/C_d$\")\n",
    "\n",
    "plt.gca().xaxis.set_major_formatter(FuncFormatter(lambda x, _: f\"{x:.1f}°\"))\n",
    "\n",
    "plt.title(\"Lift-drag ratio\")\n",
    "plt.grid()\n",
    "plt.savefig(\n",
    "    \"results/figures/01012025 - OpenFOAM - lift-drag ratio for angles of attack - random airfoil.png\",\n",
    "    dpi=300,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lift-drag vs velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_df = pd.read_csv(\"results/csv/velocity_results.csv\")\n",
    "\n",
    "velocity_df[\"velocity\"] = (\n",
    "    velocity_df[\"run_name\"].str.extract(r\"(\\d+\\.\\d+)_ms\").astype(float)\n",
    ")\n",
    "velocity_df = velocity_df[velocity_df[\"cl\"] >= 0]\n",
    "\n",
    "velocity_df = velocity_df.sort_values(by=\"velocity\")\n",
    "\n",
    "plt.plot(velocity_df[\"velocity\"], velocity_df[\"cl\"] / velocity_df[\"cd\"], \"-x\")\n",
    "plt.xlabel(\"Velocity (m/s)\")\n",
    "plt.ylabel(\"$C_l/C_d$\")\n",
    "\n",
    "plt.title(\"Lift-drag ratio\")\n",
    "plt.grid()\n",
    "plt.savefig(\n",
    "    \"results/figures/01012025 - OpenFOAM - lift-drag ratio for velocities - random airfoil.png\",\n",
    "    dpi=300,\n",
    ")\n"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}