notebooks/cylinder_drl_exercise.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "worldwide-ballet",
   "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",
    "# Controlling the flow past a cylinder with OpenFOAM and PyTorch\n",
    "\n",
    "This exercise consists of two parts, which are meant to be solved in two individual sessions. We will implement two active flow control approaches to reduce drag and lift forces acting on the surface of a circular cylinder in a stream. The base simulation is the same as in exercise 2. In both control approaches, the flow is actuated by rotating the cylinder. The rotation is easily implemented in a CFD simulation by modifying the velocity boundary condition applied to the *cylinder* patch. Say that $u_\\parallel$ is the signed cross-radial speed of the cylinder's surface. Then the surface velocity vector field $\\mathbf{u}_\\parallel$ is obtained by multiplying the speed with the tangential unit vector $\\mathbf{t}_\\parallel (\\mathbf{x})$:\n",
    "$$\n",
    "  \\mathbf{u}_\\parallel = u_\\parallel \\mathbf{t}_\\parallel (\\mathbf{x}),\n",
    "$$\n",
    "where $\\mathbf{t}_\\parallel (\\mathbf{x})$ is tangential to the cylinder's surface (note that the simulation is 2D, so there is only one sensible tangent). The two control approaches vary in the way how they determine $u_\\parallel$. In open-loop control, the control law is a function of time, e.g., $u_\\parallel = u_\\parallel (t)$. The closed-loop control law is a function of pressure sensors $\\mathbf{p}$, e.g., $u_\\parallel = u_\\parallel (\\mathbf{p})$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b45e3114",
   "metadata": {},
   "source": [
    "## Open-loop active flow control\n",
    "\n",
    "As observed in exercise 2, the uncontrolled flow past the cylinder at $Re=100$ yields almost perfectly sinusoidal lift and drag curves in the quasi-steady state (neglecting the non-zero mean). Therefore, we could try to interrupt the natural vortex shedding by rotating the cylinder periodically in alternating directions. A suitable control law might be:\n",
    "$$\n",
    "  u_\\parallel (t) = A \\mathrm{sin}(2\\pi f t),\n",
    "$$\n",
    "where the amplitude $A$ and the frequency $f$ are adjustable parameters. To implement this new velocity boundary condition, we use OpenFOAM's [inline expressions for dictionaries](https://www.openfoam.com/news/main-news/openfoam-v1912/pre-processing#pre-processing-expressions-syntax). To get started, create a dedicated folder and make a copy of the base simulation:\n",
    "```\n",
    "mkdir exercises/openloop_control\n",
    "cp -r test_cases/cylinder2D/ exercises/openloop_control/A2f6\n",
    "```\n",
    "The velocity boundary condition is specified in the *0.org/U* dictionary. Open the dictionary and make the following modifications:\n",
    "\n",
    "1. currently, there is one entry for the *cylinder*, *top*, and *bottom* patches, because the same boundary condition is applied to these patches (no-slip); remove the *cylinder* patch from the joint entry and create a new separate entry only for *cylinder*\n",
    "2. use the *exprFixedValue* type boundary condition as follows:\n",
    "```\n",
    "cylinder\n",
    "{\n",
    "  type        exprFixedValue;\n",
    "  value       $internalField;\n",
    "  amplitude   2.0;\n",
    "  frequency   6.0;\n",
    "  valueExpr   \"(vector(0, 0, 1) ^ face()) / area() * ???\";\n",
    "}\n",
    "```\n",
    "Variables can be accessed using the syntax `$variable`. E.g., `$frequency` becomes 6.0 in the example above. The expression `(vector(0, 0, 1) ^ face()) / area()` computes the tangential unit vector. The only missing piece is the expression for $u_\\parallel (t)$ in place of the three question marks. You can find the names for the required mathematical functions in the [documentation](https://www.openfoam.com/documentation/guides/latest/doc/openfoam-guide-expression-syntax.html).\n",
    "\n",
    "To run the actuated simulation, follow the usual steps:\n",
    "```\n",
    "# top-level folder of repository\n",
    "source setup-env --container\n",
    "cd exercises/openloop_control\n",
    "```\n",
    "Once the simulation is complete, open the case in ParaView:\n",
    "```\n",
    "paraview post.foam\n",
    "```\n",
    "and apply the following steps to check if the control works:\n",
    "1. set *Case Type* to *Decomposed Case*\n",
    "2. under *Mesh Regions*, select only the cylinder patch\n",
    "3. click *Apply*\n",
    "4. create a *Glyph* filter and set the *Scale Array* to *U*\n",
    "5. click the play button and observe the surface velocity field\n",
    "\n",
    "If the control is implemented correctly, we can now check its effect on the forces acting on the cylinder. To do so, create a new Jupyter notebook, load the drag and lift coefficients as in previous exercises, and plot the controlled against the un-controlled coefficients. \n",
    "\n",
    "Finally, create one or two additional parameter variations in terms of $f$ and $A$, and repeat the process of evaluating the force coefficients."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "likely-english",
   "metadata": {},
   "source": [
    "## Closed-loop active flow control\n",
    "\n",
    "The agent's aim in the DRL-based closed-loop control is learning a control law $u_\\parallel = u_\\parallel (\\mathbf{p})$ by trail and error. The state vector $\\mathbf{p}$ consists of pressure probes from 12 different locations in the cylinder's wake. To train the agent, we employ the [drlFoam](https://github.com/OFDataCommittee/drlfoam) library. This library manages the execution of various parallel simulations and the update of the agent for us. To get started, go to the drlFoam repository, follow the installation instructions, and execute the tutorial as described in the README file.\n",
    "\n",
    "The output created during the training consists of three types of files (*###* indicates the training episode):\n",
    "- *checkpoint_###.pt*: checkpoint with training details and snapshots of networks and optimizers; used to evaluate and restart a training\n",
    "- *observations_###.pt:* one or more trajectories of experience tuples based on which the agent learns\n",
    "- *policy_trace_###.pt:* snapshot of the traced policy in TorchScript format; used to evaluate a policy\n",
    "\n",
    "Meanwhile the training is running, create a Jupyter notebook that processes the *observations_###.pt* files and:\n",
    "\n",
    "1. visualize the **expected** angular velocity for various episodes\n",
    "2. visualize the drag and lift curves for various episodes\n",
    "3. visualize the mean reward over the episodes\n",
    "\n",
    "The *observations* are stored as simple Pickle files and can be loaded with `torch.load(\"file_name.pt\")`. Each file contains a list of trajectories. Individual trajectories are organized as dictionaries. The following lines of code demonstrate how to assess the available trajectory content:\n",
    "\n",
    "```\n",
    "import torch as pt\n",
    "# zeroth episode\n",
    "obs = pt.load(\"observations_0.pt\")\n",
    "# second trajectory\n",
    "tr_1 = obs[1]\n",
    "# trajectory content\n",
    "print(tr_1.keys())\n",
    "```\n",
    "\n",
    "Once the training is complete, evaluate the best policy according to the mean rewards:\n",
    "- create a copy of the *base* simulation (inside the training folder)\n",
    "- copy the best policy trace into the copy of *base* and rename it to *policy.pt*\n",
    "- in the *U* file, set *train* to *false*\n",
    "- in the *controlDict*  \n",
    "  - set the *endTime* to 8  \n",
    "  - set *purgeWrite* to 0  \n",
    "  - set *writeInterval* to 0.01  \n",
    "  - set *timeStart* to 0 for both *forces* and *probes* function object\n",
    "- execute the simulation with *Allrun.pre* and *Allrun*\n",
    "- once the simulation is complete, inspect the results in ParaView\n",
    "\n",
    "**Congratulations! This completes the final two exercise session.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7525d54",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 ('ml-cfd': venv)",
   "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "cbb112329e189d437c2dd20cb32069784312f3fb0cc6892ecadf8f4f78f994e7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}