Making model parameters in demos consistent and setting default seed.

Author: matt-graham

pompy/demos.py

@@ -15,7 +15,10 @@
 from pompy import models, processors
 
 
-def set_up_figure(fig_size=(8, 4)):
+DEFAULT_SEED = 20181108
+
+
+def set_up_figure(fig_size=(10, 5)):
     """Set up Matplotlib figure with simulation time title text.
 
     Parameters
@@ -44,7 +47,7 @@ def wrapped_update(i):
     return inner_decorator
 
 
-def wind_model_demo(dt=0.01, t_max=100, steps_per_frame=20):
+def wind_model_demo(dt=0.01, t_max=100, steps_per_frame=20, seed=DEFAULT_SEED):
     """Set up wind model and animate velocity field with quiver plot.
 
     Parameters
@@ -55,6 +58,8 @@ def wind_model_demo(dt=0.01, t_max=100, steps_per_frame=20):
         End time to simulate to.
     steps_per_frame: integer
         Number of simulation time steps to perform between animation frames.
+    seed : integer
+        Seed for random number generator.
 
     Returns
     -------
@@ -65,10 +70,14 @@ def wind_model_demo(dt=0.01, t_max=100, steps_per_frame=20):
     anim : FuncAnimation
         Matplotlib animation object.
     """
+    rng = np.random.RandomState(seed)
     # define simulation region
-    wind_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)
+    wind_region = models.Rectangle(x_min=0., x_max=100., y_min=-25., y_max=25.)
     # set up wind model
-    wind_model = models.WindModel(wind_region, 21, 11, u_av=1., k_x=2., k_y=2.)
+    wind_model = models.WindModel(wind_region, 21, 11, rng=rng)
+    # let simulation run for 10s to equilibrate wind model
+    for t in np.arange(0, 10, dt):
+        wind_model.update(dt)
     # generate figure and attach close event
     fig, ax, title = set_up_figure()
     # create quiver plot of initial velocity field
@@ -95,7 +104,8 @@ def update(i):
     return fig, ax, anim
 
 
-def plume_model_demo(dt=0.01, t_max=100, steps_per_frame=200):
+def plume_model_demo(dt=0.01, t_max=100, steps_per_frame=200,
+                     seed=DEFAULT_SEED):
     """Set up plume model and animate puffs overlayed over velocity field.
 
     Puff positions displayed using Matplotlib `scatter` plot function and
@@ -110,6 +120,8 @@ def plume_model_demo(dt=0.01, t_max=100, steps_per_frame=200):
         End time to simulate to.
     steps_per_frame: integer
         Number of simulation time steps to perform between animation frames.
+    seed : integer
+        Seed for random number generator.
 
     Returns
     -------
@@ -120,15 +132,17 @@ def plume_model_demo(dt=0.01, t_max=100, steps_per_frame=200):
     anim : FuncAnimation
         Matplotlib animation object.
     """
+    rng = np.random.RandomState(seed)
     # define simulation region
-    wind_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)
+    sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)
     # set up wind model
-    wind_model = models.WindModel(wind_region, 21, 11, u_av=1., k_x=2., k_y=2.)
+    wind_model = models.WindModel(sim_region, 21, 11, rng=rng)
+    # let simulation run for 10s to equilibrate wind model
+    for t in np.arange(0, 10, dt):
+        wind_model.update(dt)
     # set up plume model
-    sim_region = models.Rectangle(x_min=0., x_max=50, y_min=-12.5, y_max=12.5)
     plume_model = models.PlumeModel(
-        sim_region, (5., 0., 0.), wind_model, puff_release_rate=10,
-        centre_rel_diff_scale=2)
+        sim_region, (5., 0., 0.), wind_model, rng=rng)
     # set up figure window
     fig, ax, title = set_up_figure()
     # create quiver plot of initial velocity field
@@ -168,12 +182,14 @@ def update(i):
     return fig, ax, anim
 
 
-def concentration_array_demo(dt=0.01, t_max=100, steps_per_frame=50):
-    """Set up plume model and animate concentration fields.
+def conc_point_val_demo(dt=0.01, t_max=5, steps_per_frame=1, x=10., y=0.0,
+                        seed=DEFAULT_SEED):
+    """Set up plume model and animate concentration at a point as time series.
 
     Demonstration of setting up plume model and processing the outputted
-    puff arrays with the `ConcentrationArrayGenerator` class, the resulting
-    arrays being displayed with the Matplotlib `imshow` function.
+    puff arrays with the ConcentrationPointValueCalculator class, the
+    resulting concentration time course at a point in the odour plume being
+    displayed with the Matplotlib `plot` function.
 
     Parameters
     ----------
@@ -183,6 +199,12 @@ def concentration_array_demo(dt=0.01, t_max=100, steps_per_frame=50):
         End time to simulate to.
     steps_per_frame: integer
         Number of simulation time steps to perform between animation frames.
+    x : float
+        x-coordinate of point to measure concentration at.
+    y : float
+        y-coordinate of point to measure concentration at.
+    seed : integer
+        Seed for random number generator.
 
     Returns
     -------
@@ -193,49 +215,56 @@ def concentration_array_demo(dt=0.01, t_max=100, steps_per_frame=50):
     anim : FuncAnimation
         Matplotlib animation object.
     """
+    rng = np.random.RandomState(seed)
     # define simulation region
-    wind_region = models.Rectangle(x_min=0., x_max=10., y_min=-2., y_max=2.)
-    sim_region = models.Rectangle(x_min=0., x_max=4, y_min=-1., y_max=1.)
+    sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)
     # set up wind model
-    wind_model = models.WindModel(wind_region, 21, 11, u_av=1.,)
+    wind_model = models.WindModel(sim_region, 21, 11, rng=rng)
     # set up plume model
     plume_model = models.PlumeModel(
-        sim_region, (0.1, 0., 0.), wind_model, centre_rel_diff_scale=1.5,
-        puff_release_rate=500, puff_init_rad=0.001)
-    # set up concentration array generator
-    array_gen = processors.ConcentrationArrayGenerator(
-        sim_region, 0.01, 500, 500, 1.)
+        sim_region, (5., 0., 0.), wind_model, rng=rng)
+    # let simulation run for 10s to initialise models
+    for t in np.arange(0, 10, dt):
+        wind_model.update(dt)
+        plume_model.update(dt)
+    # set up concentration point value calculator
+    val_calc = processors.ConcentrationValueCalculator(1.)
+    conc_vals = []
+    conc_vals.append(val_calc.calc_conc_point(plume_model.puff_array, x, y))
+    ts = [0.]
     # set up figure
     fig, ax, title = set_up_figure()
     # display initial concentration field as image
-    conc_array = array_gen.generate_single_array(plume_model.puff_array)
-    conc_im = plt.imshow(conc_array.T, extent=sim_region, vmin=0, vmax=3e4,
-                         cmap='Reds')
-    ax.set_xlabel('x-coordinate / m')
-    ax.set_ylabel('y-coordinate / m')
-    ax.set_aspect(1)
+    conc_line, = plt.plot(ts, conc_vals)
+    ax.set_xlim(0., t_max)
+    ax.set_ylim(0., 150.)
+    ax.set_xlabel('Time / s')
+    ax.set_ylabel('Normalised concentration')
+    ax.grid(True)
     fig.tight_layout()
 
     # define update function
     @update_decorator(dt, title, steps_per_frame, [wind_model, plume_model])
     def update(i):
-        conc_im.set_data(
-            array_gen.generate_single_array(plume_model.puff_array).T)
-        return [conc_im]
+        ts.append(dt * i * steps_per_frame)
+        conc_vals.append(
+            val_calc.calc_conc_point(plume_model.puff_array, x, y))
+        conc_line.set_data(ts, conc_vals)
+        return [conc_line]
 
     # create animation object
     n_frame = int(t_max / (dt * steps_per_frame) + 0.5)
     anim = FuncAnimation(fig, update, frames=n_frame, blit=True)
     return fig, ax, anim
 
 
-def conc_point_val_demo(dt=0.01, t_max=5, steps_per_frame=1, x=1., y=0.0):
-    """Set up plume model and animate concentration at a point as time series.
+def concentration_array_demo(dt=0.01, t_max=100, steps_per_frame=50,
+                             seed=DEFAULT_SEED):
+    """Set up plume model and animate concentration fields.
 
     Demonstration of setting up plume model and processing the outputted
-    puff arrays with the ConcentrationPointValueCalculator class, the
-    resulting concentration time course at a point in the odour plume being
-    displayed with the Matplotlib `plot` function.
+    puff arrays with the `ConcentrationArrayGenerator` class, the resulting
+    arrays being displayed with the Matplotlib `imshow` function.
 
     Parameters
     ----------
@@ -245,10 +274,8 @@ def conc_point_val_demo(dt=0.01, t_max=5, steps_per_frame=1, x=1., y=0.0):
         End time to simulate to.
     steps_per_frame: integer
         Number of simulation time steps to perform between animation frames.
-    x : float
-        x-coordinate of point to measure concentration at.
-    y : float
-        y-coordinate of point to measure concentration at.
+    seed : integer
+        Seed for random number generator.
 
     Returns
     -------
@@ -259,43 +286,38 @@ def conc_point_val_demo(dt=0.01, t_max=5, steps_per_frame=1, x=1., y=0.0):
     anim : FuncAnimation
         Matplotlib animation object.
     """
+    rng = np.random.RandomState(seed)
     # define simulation region
-    wind_region = models.Rectangle(x_min=0., x_max=10., y_min=-2., y_max=2.)
-    sim_region = models.Rectangle(x_min=0., x_max=2, y_min=-1., y_max=1.)
+    sim_region = models.Rectangle(x_min=0., x_max=100, y_min=-25., y_max=25.)
     # set up wind model
-    wind_model = models.WindModel(wind_region, 21, 11, u_av=2.)
+    wind_model = models.WindModel(sim_region, 21, 11, rng=rng)
     # set up plume model
     plume_model = models.PlumeModel(
-        sim_region, (0.1, 0., 0.), wind_model, centre_rel_diff_scale=1.5,
-        puff_release_rate=25, puff_init_rad=0.01)
-    # let simulation run for 10s to get plume established
+        sim_region, (5., 0., 0.), wind_model, rng=rng)
+    # let simulation run for 10s to initialise models
     for t in np.arange(0, 10, dt):
         wind_model.update(dt)
         plume_model.update(dt)
-    # set up concentration point value calculator
-    val_calc = processors.ConcentrationValueCalculator(0.01**2)
-    conc_vals = []
-    conc_vals.append(val_calc.calc_conc_point(plume_model.puff_array, x, y))
-    ts = [0.]
+    # set up concentration array generator
+    array_gen = processors.ConcentrationArrayGenerator(
+        sim_region, 0.01, 500, 250, 1.)
     # set up figure
     fig, ax, title = set_up_figure()
     # display initial concentration field as image
-    conc_line, = plt.plot(ts, conc_vals)
-    ax.set_xlim(0., t_max)
-    ax.set_ylim(0., .5)
-    ax.set_xlabel('Time / s')
-    ax.set_ylabel('Normalised concentration')
-    ax.grid(True)
+    conc_array = array_gen.generate_single_array(plume_model.puff_array)
+    conc_im = plt.imshow(conc_array.T, extent=sim_region, cmap='Reds',
+                         vmin=0., vmax=1.)
+    ax.set_xlabel('x-coordinate / m')
+    ax.set_ylabel('y-coordinate / m')
+    ax.set_aspect(1)
     fig.tight_layout()
 
     # define update function
     @update_decorator(dt, title, steps_per_frame, [wind_model, plume_model])
     def update(i):
-        ts.append(dt * i * steps_per_frame)
-        conc_vals.append(
-            val_calc.calc_conc_point(plume_model.puff_array, x, y))
-        conc_line.set_data(ts, conc_vals)
-        return [conc_line]
+        conc_im.set_data(
+            array_gen.generate_single_array(plume_model.puff_array).T)
+        return [conc_im]
 
     # create animation object
     n_frame = int(t_max / (dt * steps_per_frame) + 0.5)