Simple example of dynamic control

Import modules

import numpy as np
import pandas as pd
from pipedream_solver.hydraulics import SuperLink
from pipedream_solver.simulation import Simulation

import matplotlib.pyplot as plt
import seaborn as sns

Load model data

input_path = '../data/hillslope'
superjunctions = pd.read_csv(f'{input_path}/hillslope_control_superjunctions.csv')
superlinks = pd.read_csv(f'{input_path}/hillslope_control_superlinks.csv')

superjunctions

	name	id	z_inv	h_0	bc	storage	a	b	c	max_depth	map_x	map_y
0	0	0	2	0.00001	False	functional	0.0	0.0	200.0	inf	0	0
1	1	1	1	0.00001	False	functional	0.0	0.0	1000.0	inf	1	1
2	2	2	1	0.00001	False	functional	0.0	0.0	1000.0	inf	2	2
3	3	3	0	0.00001	False	functional	0.0	0.0	1000.0	inf	3	3

superlinks

	name	id	sj_0	sj_1	in_offset	out_offset	dx	n	shape	g1	g2	g3	g4	Q_0	h_0	ctrl	A_s	A_c	C
0	0	0	0	1	0.0	0.0	1000	0.035	rect_open	10	5	0	0	0	0.00001	False	100	0	0
1	0	0	2	3	0.0	0.0	1000	0.035	rect_open	10	5	0	0	0	0.00001	False	100	0	0

Create a control orifice

orifices = {
    0 : {
        'id' : 0,
        'sj_0' : 1,
        'sj_1' : 2,
        'A' : 0.3048**2,
        'orientation' : 'side',
        'z_o' : 0,
        'y_max' : 0.3048,
        'C' : 0.67
    }
}

orifices = pd.DataFrame.from_dict(orifices, orient='index')

orifices

	id	sj_0	sj_1	A	orientation	z_o	y_max	C
0	0	1	2	0.092903	side	0	0.3048	0.67

Instantiate model

internal_links = 24
model = SuperLink(superlinks, superjunctions, orifices=orifices,
                  internal_links=internal_links)
model.spinup(n_steps=1000)

dt = 10                            # Model time step (s)
Q_in = 1e-3 * np.ones(model.M)     # Flow into each internal junction (cms)
Q_0Ik = 1e-3 * np.ones(model.NIk)  # Flow into each internal junction (cms)
u_0 = np.zeros(len(orifices))      # Orifice closed signal
u_1 = np.ones(len(orifices))       # Orifice open signal

# Plot profile from superjunctions 0 to 1 (uphill to downhill)
_ = model.plot_profile([0, 1, 2, 3], width=100)

Run model with orifice closed

# For each timestep...
while model.t < (8 * 3600):
    if model.t < (4 * 3600):
        model.step(dt=dt, Q_in=Q_in, Q_0Ik=Q_0Ik, u_o=u_0)
    else:
        model.step(dt=dt, u_o=u_0)

_ = model.plot_profile([0, 1, 2, 3], width=100)

Continue running model with orifice open

while model.t < (12 * 3600):
    model.step(dt=dt, u_o=u_1)

_ = model.plot_profile([0, 1, 2, 3], width=100)

Run full simulation with control

model = SuperLink(superlinks, superjunctions, orifices=orifices,
                  internal_links=internal_links)
model.spinup(n_steps=1000)

# Create simulation context manager
with Simulation(model, t_start=0, t_end=(24 * 3600)) as simulation:
    # While simulation time has not expired...
    while simulation.t <= simulation.t_end:
        if simulation.t < (8 * 3600):
            simulation.step(dt=dt, Q_in=Q_in, Q_Ik=Q_0Ik, u_o=u_0)
        elif (simulation.t > 8 * 3600) and (simulation.t < 12 * 3600):
            simulation.step(dt=dt, u_o=u_0)
        else:
            simulation.step(dt=dt, u_o=u_1)
        # Record internal depth and flow states
        simulation.record_state()
        # Reposition junctions
        simulation.model.reposition_junctions()
        # Print progress bar
        simulation.print_progress()

[==================================================] 100.0% [4.77 s]

Visualize results

Profile at end of simulation

_ = model.plot_profile([0, 1, 2, 3], width=100)

Superjunction heads

sns.set_palette('husl')
simulation.states.H_j.plot()
plt.axvline(12 * 3600, c='k', alpha=0.3, linestyle='--')
plt.title('Superjunction heads')
plt.ylabel('Head (m)')
_ = plt.xlabel('Time (s)')

Internal superlink depths

fig, ax = plt.subplots(2, figsize=(10, 10))

(simulation.states.h_Ik[np.flatnonzero(model._kI == 0)]
 .plot(ax=ax[0], legend=False, alpha=0.8, cmap='rainbow_r'))
(simulation.states.h_Ik[np.flatnonzero(model._kI == 1)]
 .plot(ax=ax[1], legend=False, alpha=0.8, cmap='rainbow_r'))

ax[0].axvline(12 * 3600, c='k', alpha=0.3, linestyle='--')
ax[1].axvline(12 * 3600, c='k', alpha=0.3, linestyle='--')

ax[0].set_title('Depth in internal junctions (upstream of orifice)')
ax[1].set_title('Depth in internal junctions (downstream of orifice)')
ax[0].xaxis.set_ticklabels([])
ax[0].set_ylabel('Depth (m)')
ax[1].set_ylabel('Depth (m)')
ax[1].set_xlabel('Time (s)')

Internal superlink flows

fig, ax = plt.subplots(2, figsize=(10, 10))

(simulation.states.Q_ik[np.flatnonzero(model._ki == 0)]
 .plot(ax=ax[0], legend=False, alpha=0.8, cmap='rainbow_r'))
(simulation.states.Q_ik[np.flatnonzero(model._ki == 1)]
 .plot(ax=ax[1], legend=False, alpha=0.8, cmap='rainbow_r'))

ax[0].axvline(12 * 3600, c='k', alpha=0.3, linestyle='--')
ax[1].axvline(12 * 3600, c='k', alpha=0.3, linestyle='--')

ax[0].set_title('Flow in internal links (upstream of orifice)')
ax[1].set_title('Flow in internal links (downstream of orifice)')
ax[0].xaxis.set_ticklabels([])
ax[0].set_ylabel('Flow (cms)')
ax[1].set_ylabel('Flow (cms)')
ax[1].set_xlabel('Time (s)')