{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from pompy import demos\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams['animation.html'] = 'html5'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# `pompy` demos\n",
"\n",
"A number of demonstrations of setting up plume models with `pompy` are included in the `pompy.demos` modules. This notebook displays the code of these demonstration function along side cells allowing the demos to be run and the resulting visualisations displayed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting wind model velocity field\n",
"\n",
"Simulate a wind velocity field through time plotting the velocity fields as an animated quiver plot.\n",
"\n",
"```Python\n",
"def wind_model_demo(dt=0.01, t_max=100, steps_per_frame=20, seed=DEFAULT_SEED):\n",
" \"\"\"Set up wind model and animate velocity field with quiver plot.\n",
"\n",
" Parameters\n",
" ----------\n",
" dt : float\n",
" Simulation timestep.\n",
" t_max : float\n",
" End time to simulate to.\n",
" steps_per_frame: integer\n",
" Number of simulation time steps to perform between animation frames.\n",
" seed : integer\n",
" Seed for random number generator.\n",
"\n",
" Returns\n",
" -------\n",
" fig : Figure\n",
" Matplotlib figure object.\n",
" ax : AxesSubplot\n",
" Matplotlib axis object.\n",
" anim : FuncAnimation\n",
" Matplotlib animation object.\n",
" \"\"\"\n",
" rng = np.random.RandomState(seed)\n",
" # define simulation region\n",
" wind_region = models.Rectangle(x_min=0., x_max=100., y_min=-25., y_max=25.)\n",
" # set up wind model\n",
" wind_model = models.WindModel(wind_region, 21, 11, rng=rng)\n",
" # let simulation run for 10s to equilibrate wind model\n",
" for t in np.arange(0, 10, dt):\n",
" wind_model.update(dt)\n",
" # generate figure and attach close event\n",
" fig, ax, title = set_up_figure()\n",
" # create quiver plot of initial velocity field\n",
" vf_plot = ax.quiver(wind_model.x_points, wind_model.y_points,\n",
" wind_model.velocity_field.T[0],\n",
" wind_model.velocity_field.T[1], width=0.003)\n",
" # expand axis limits to make vectors at boundary of field visible\n",
" ax.axis(ax.axis() + np.array([-0.25, 0.25, -0.25, 0.25]))\n",
" ax.set_xlabel('x-coordinate / m')\n",
" ax.set_ylabel('y-coordinate / m')\n",
" ax.set_aspect(1)\n",
" fig.tight_layout()\n",
"\n",
" # define update function\n",
" @update_decorator(dt, title, steps_per_frame, [wind_model])\n",
" def update(i):\n",
" vf_plot.set_UVC(\n",
" wind_model.velocity_field.T[0], wind_model.velocity_field.T[1])\n",
" return [vf_plot]\n",
"\n",
" # create animation object\n",
" n_frame = int(t_max / (dt * steps_per_frame) + 0.5)\n",
" anim = FuncAnimation(fig, update, n_frame, blit=True)\n",
" return fig, ax, anim\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax, anim = demos.wind_model_demo()\n",
"plt.close(fig)\n",
"anim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Combining wind and plume model visualisation\n",
"\n",
"Simulate both a wind and odour puff plume model, visualising puffs using an animated scatter plot overlayed over a quiver plot of wind velocity field.\n",
"\n",
"\n",
"```Python\n",
"def plume_model_demo(dt=0.01, t_max=100, steps_per_frame=200,\n",
" seed=DEFAULT_SEED):\n",
" \"\"\"Set up plume model and animate puffs overlayed over velocity field.\n",
"\n",
" Puff positions displayed using Matplotlib `scatter` plot function and\n",
" velocity field displayed using `quiver` plot function.\n",
" plot and quiver functions.\n",
"\n",
" Parameters\n",
" ----------\n",
" dt : float\n",
" Simulation timestep.\n",
" t_max : float\n",
" End time to simulate to.\n",
" steps_per_frame: integer\n",
" Number of simulation time steps to perform between animation frames.\n",
" seed : integer\n",
" Seed for random number generator.\n",
"\n",
" Returns\n",
" -------\n",
" fig : Figure\n",
" Matplotlib figure object.\n",
" ax : AxesSubplot\n",
" Matplotlib axis object.\n",
" anim : FuncAnimation\n",
" Matplotlib animation object.\n",
" \"\"\"\n",
" rng = np.random.RandomState(seed)\n",
" # define simulation region\n",
" sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)\n",
" # set up wind model\n",
" wind_model = models.WindModel(sim_region, 21, 11, rng=rng)\n",
" # let simulation run for 10s to equilibrate wind model\n",
" for t in np.arange(0, 10, dt):\n",
" wind_model.update(dt)\n",
" # set up plume model\n",
" plume_model = models.PlumeModel(\n",
" sim_region, (5., 0., 0.), wind_model, rng=rng)\n",
" # set up figure window\n",
" fig, ax, title = set_up_figure()\n",
" # create quiver plot of initial velocity field\n",
" # quiver expects first array dimension (rows) to correspond to y-axis\n",
" # therefore need to transpose\n",
" vf_plot = plt.quiver(\n",
" wind_model.x_points, wind_model.y_points,\n",
" wind_model.velocity_field.T[0], wind_model.velocity_field.T[1],\n",
" width=0.003)\n",
" # expand axis limits to make vectors at boundary of field visible\n",
" ax.axis(ax.axis() + np.array([-0.25, 0.25, -0.25, 0.25]))\n",
" # draw initial puff positions with scatter plot\n",
" radius_mult = 200\n",
" pp_plot = plt.scatter(\n",
" plume_model.puff_array[:, 0], plume_model.puff_array[:, 1],\n",
" radius_mult * plume_model.puff_array[:, 3]**0.5, c='r',\n",
" edgecolors='none')\n",
" ax.set_xlabel('x-coordinate / m')\n",
" ax.set_ylabel('y-coordinate / m')\n",
" ax.set_aspect(1)\n",
" fig.tight_layout()\n",
"\n",
" # define update function\n",
" @update_decorator(dt, title, steps_per_frame, [wind_model, plume_model])\n",
" def update(i):\n",
" # update velocity field quiver plot data\n",
" vf_plot.set_UVC(wind_model.velocity_field[:, :, 0].T,\n",
" wind_model.velocity_field[:, :, 1].T)\n",
" # update puff position scatter plot positions and sizes\n",
" pp_plot.set_offsets(plume_model.puff_array[:, :2])\n",
" pp_plot._sizes = radius_mult * plume_model.puff_array[:, 3]**0.5\n",
" return [vf_plot, pp_plot]\n",
"\n",
" # create animation object\n",
" n_frame = int(t_max / (dt * steps_per_frame) + 0.5)\n",
" anim = FuncAnimation(fig, update, frames=n_frame, blit=True)\n",
" return fig, ax, anim\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax, anim = demos.plume_model_demo()\n",
"plt.close(fig)\n",
"anim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting the variation in concentration at a point\n",
"\n",
"Simulate a plume model and compute concentration at a point in simulation region, plotting as an animate time series.\n",
"\n",
"```Python\n",
"def conc_point_val_demo(dt=0.01, t_max=5, steps_per_frame=1, x=10., y=0.0,\n",
" seed=DEFAULT_SEED):\n",
" \"\"\"Set up plume model and animate concentration at a point as time series.\n",
"\n",
" Demonstration of setting up plume model and processing the outputted\n",
" puff arrays with the ConcentrationPointValueCalculator class, the\n",
" resulting concentration time course at a point in the odour plume being\n",
" displayed with the Matplotlib `plot` function.\n",
"\n",
" Parameters\n",
" ----------\n",
" dt : float\n",
" Simulation timestep.\n",
" t_max : float\n",
" End time to simulate to.\n",
" steps_per_frame: integer\n",
" Number of simulation time steps to perform between animation frames.\n",
" x : float\n",
" x-coordinate of point to measure concentration at.\n",
" y : float\n",
" y-coordinate of point to measure concentration at.\n",
" seed : integer\n",
" Seed for random number generator.\n",
"\n",
" Returns\n",
" -------\n",
" fig : Figure\n",
" Matplotlib figure object.\n",
" ax : AxesSubplot\n",
" Matplotlib axis object.\n",
" anim : FuncAnimation\n",
" Matplotlib animation object.\n",
" \"\"\"\n",
" rng = np.random.RandomState(seed)\n",
" # define simulation region\n",
" sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)\n",
" # set up wind model\n",
" wind_model = models.WindModel(sim_region, 21, 11, rng=rng)\n",
" # set up plume model\n",
" plume_model = models.PlumeModel(\n",
" sim_region, (5., 0., 0.), wind_model, rng=rng)\n",
" # let simulation run for 10s to initialise models\n",
" for t in np.arange(0, 10, dt):\n",
" wind_model.update(dt)\n",
" plume_model.update(dt)\n",
" # set up concentration point value calculator\n",
" val_calc = processors.ConcentrationValueCalculator(1.)\n",
" conc_vals = []\n",
" conc_vals.append(val_calc.calc_conc_point(plume_model.puff_array, x, y))\n",
" ts = [0.]\n",
" # set up figure\n",
" fig, ax, title = set_up_figure()\n",
" # display initial concentration field as image\n",
" conc_line, = plt.plot(ts, conc_vals)\n",
" ax.set_xlim(0., t_max)\n",
" ax.set_ylim(0., 150.)\n",
" ax.set_xlabel('Time / s')\n",
" ax.set_ylabel('Normalised concentration')\n",
" ax.grid(True)\n",
" fig.tight_layout()\n",
"\n",
" # define update function\n",
" @update_decorator(dt, title, steps_per_frame, [wind_model, plume_model])\n",
" def update(i):\n",
" ts.append(dt * i * steps_per_frame)\n",
" conc_vals.append(\n",
" val_calc.calc_conc_point(plume_model.puff_array, x, y))\n",
" conc_line.set_data(ts, conc_vals)\n",
" return [conc_line]\n",
"\n",
" # create animation object\n",
" n_frame = int(t_max / (dt * steps_per_frame) + 0.5)\n",
" anim = FuncAnimation(fig, update, frames=n_frame, blit=True)\n",
" return fig, ax, anim\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax, anim = demos.conc_point_val_demo()\n",
"plt.close(fig)\n",
"anim"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualising concentration fields"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simulate a plume model and visualise odour concentration fields as animated images.\n",
"\n",
"```Python\n",
"def concentration_array_demo(dt=0.01, t_max=100, steps_per_frame=50,\n",
" seed=DEFAULT_SEED):\n",
" \"\"\"Set up plume model and animate concentration fields.\n",
"\n",
" Demonstration of setting up plume model and processing the outputted\n",
" puff arrays with the `ConcentrationArrayGenerator` class, the resulting\n",
" arrays being displayed with the Matplotlib `imshow` function.\n",
"\n",
" Parameters\n",
" ----------\n",
" dt : float\n",
" Simulation timestep.\n",
" t_max : float\n",
" End time to simulate to.\n",
" steps_per_frame: integer\n",
" Number of simulation time steps to perform between animation frames.\n",
" seed : integer\n",
" Seed for random number generator.\n",
"\n",
" Returns\n",
" -------\n",
" fig : Figure\n",
" Matplotlib figure object.\n",
" ax : AxesSubplot\n",
" Matplotlib axis object.\n",
" anim : FuncAnimation\n",
" Matplotlib animation object.\n",
" \"\"\"\n",
" rng = np.random.RandomState(seed)\n",
" # define simulation region\n",
" sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)\n",
" # set up wind model\n",
" wind_model = models.WindModel(sim_region, 21, 11, rng=rng)\n",
" # set up plume model\n",
" plume_model = models.PlumeModel(\n",
" sim_region, (5., 0., 0.), wind_model, rng=rng)\n",
" # let simulation run for 10s to initialise models\n",
" for t in np.arange(0, 10, dt):\n",
" wind_model.update(dt)\n",
" plume_model.update(dt)\n",
" # set up concentration array generator\n",
" array_gen = processors.ConcentrationArrayGenerator(\n",
" sim_region, 0.01, 500, 250, 1.)\n",
" # set up figure\n",
" fig, ax, title = set_up_figure()\n",
" # display initial concentration field as image\n",
" conc_array = array_gen.generate_single_array(plume_model.puff_array)\n",
" conc_im = plt.imshow(conc_array.T, extent=sim_region, cmap='Reds',\n",
" vmin=0., vmax=1.)\n",
" ax.set_xlabel('x-coordinate / m')\n",
" ax.set_ylabel('y-coordinate / m')\n",
" ax.set_aspect(1)\n",
" fig.tight_layout()\n",
"\n",
" # define update function\n",
" @update_decorator(dt, title, steps_per_frame, [wind_model, plume_model])\n",
" def update(i):\n",
" conc_im.set_data(\n",
" array_gen.generate_single_array(plume_model.puff_array).T)\n",
" return [conc_im]\n",
"\n",
" # create animation object\n",
" n_frame = int(t_max / (dt * steps_per_frame) + 0.5)\n",
" anim = FuncAnimation(fig, update, frames=n_frame, blit=True)\n",
" return fig, ax, anim\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, ax, anim = demos.concentration_array_demo()\n",
"plt.close(fig)\n",
"anim"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}