Commit d1bc9d42 authored by Jesse VanderWees's avatar Jesse VanderWees 🐘 Committed by Jorn Baayen

Add an example demonstrating the different flavours of the Saint-Venant...

Add an example demonstrating the different flavours of the Saint-Venant equation implemented in RTC-Tools. The results are compared with results from HEC-RAS.
parent 9b3d5df2
......@@ -370,6 +370,7 @@ MODELICA_EXAMPLES = [
"lookup_table",
"simulation",
"cascading_channels",
"channel_pulse",
"channel_wave_damping",
]
......
......@@ -12,4 +12,5 @@ This section provides examples demonstrating key features of RTC-Tools optimizat
optimization/lookup_table
optimization/ensemble
optimization/cascading_channels
optimization/channel_pulse
optimization/channel_wave_damping
Modeling Waves in Rivers and Canals
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. image:: ../../images/Lake-Hartwell-24012016.jpg
.. :href: https://pixabay.com/en/pont-du-gard-aqueduct-roman-unesco-1742029/
.. pixabay content is released under a CC0 Public Domain licence - no attribution needed
.. note::
This is a more advanced example that implements advanced channel hydraulics
in RTC-Tools. It also capitalizes on the homotopy techniques available in
RTC-Tools. If you are a first-time user of RTC-Tools, see :doc:`basic`.
The RTC-Tools is capable of handling non-linear hydraulics. In this example, we
model a river channel that receives a sudden pulse of higher-than-usual water
volumes. We compare the results to those of an identical model built in HEC-RAS.
The Model
---------
In this example, water is flowing through a single channel. There is an inflow
at the upstream end and a water level bound at the downstream end.
In OpenModelica Connection Editor, the model looks like this in plan view:
.. image:: ../../images/single_channel.png
In text mode, the Modelica model looks as follows (with annotation statements
removed):
.. literalinclude:: ../../_build/mo/channel_pulse.mo
:language: modelica
:lineno-match:
The plan view of the model looks like this in HEC-RAS:
.. image:: ../../images/single_channel_hec-ras.png
The channel cross-section is a simple trapezoidal shape. As rendered by HEC-RAS,
here is a cross-section view of the channel being modeled:
.. image:: ../../images/channel_impulse_xs.png
The model was built with HEC-RAS version 5.0.6. In case you wish to verify the
HEC-RAS model yourself, a zip of the HEC-RAS model used in this comparison is
available: :download:`HEC-RAS.zip <../../../examples/channel_pulse/HEC-RAS.zip>`
The Python File
---------------
To keep this example simple and to allow for a 1:1 comparison with HEC-RAS, we
will not have any decision variables in this model.
.. literalinclude:: ../../../examples/channel_pulse/src/example.py
:language: python
:lineno-match:
As you can see, this model is as simple as it gets. We only add a constraint to
keep the initialization states consistent with the HEC-RAS initialization.
Comparison of Numerical Schemes
-------------------------------
HEC-RAS and RTC-Tools use different numerical schemes:
* RTC-Tools uses a staggered grid, whereas HEC-RAS uses a box scheme;
* RTC-Tools provides a solution of the original nonlinear equations, whereas
HEC-RAS uses a linearization technique to compute discharge and level
increments.
Comparison of Results
---------------------
The results from the RTC-Tools run are found in the output directory with the
name ``timeseries_export.csv``, and the results generated by HEC-RAS have been
exported into the same directory under the name ``HEC-RAS_results.csv``. We can
compare the results using the Python library ``matplotlib``:
.. plot:: examples/pyplots/channel_pulse_results.py
from datetime import datetime
from pathlib import Path
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
# Map output_dir
output_dir = Path("../../../examples/channel_pulse/output/").resolve()
# Import Data
rtc_tools_record = np.recfromcsv(output_dir / "timeseries_export_inertial_wave.csv", encoding=None)
rtc_tools_conv_acc_record = np.recfromcsv(
output_dir / "timeseries_export_saint_venant.csv", encoding=None
)
rtc_tools_conv_acc_upwind_record = np.recfromcsv(
output_dir / "timeseries_export_saint_venant_upwind.csv", encoding=None
)
hec_ras_record = np.recfromcsv(output_dir / "HEC-RAS_export.csv", encoding=None)
# Get times as datetime objects
times = list(
map(lambda x: datetime.strptime(x, "%Y-%m-%d %H:%M:%S"), rtc_tools_record["time"])
)
# Generate Plot
n_subplots = 2
fig, axarr = plt.subplots(n_subplots, sharex=True, figsize=(8, 4 * n_subplots))
axarr[0].set_title("Water Levels and Flow Rates")
# Upper subplot
axarr[0].set_ylabel("Flow Rate [m³/s]")
axarr[0].plot(
times,
rtc_tools_record["channel_q_up"],
label="Upstream",
color="xkcd:dark sky blue",
)
axarr[0].plot(
times,
rtc_tools_record["channel_q_dn"],
label="Downstream\n(RTC-Tools Inertial Wave)",
linestyle="--",
color="red",
)
axarr[0].plot(
times,
rtc_tools_conv_acc_record["channel_q_dn"],
label="Downstream\n(RTC-Tools Saint Venant central diff.)",
linestyle="--",
color="darkorange",
)
axarr[0].plot(
times,
rtc_tools_conv_acc_upwind_record["channel_q_dn"],
label="Downstream\n(RTC-Tools Saint Venant upwind)",
linestyle="--",
color="purple",
)
axarr[0].plot(
times,
hec_ras_record["channel_q_dn"],
label="Downstream\n(HEC-RAS)",
linestyle="--",
color="darkgreen",
)
# Lower subplot
axarr[1].set_ylabel("Water Level [m]")
axarr[1].plot(
times,
rtc_tools_record["channel_h_up"],
label="Upstream\n(RTC-Tools Inertial Wave)",
linestyle="--",
color="red",
)
axarr[1].plot(
times,
rtc_tools_conv_acc_record["channel_h_up"],
label="Upstream\n(RTC-Tools Saint Venant central diff.)",
linestyle="--",
color="darkorange",
)
axarr[1].plot(
times,
rtc_tools_conv_acc_upwind_record["channel_h_up"],
label="Upstream\n(RTC-Tools Saint Venant upwind)",
linestyle="--",
color="purple",
)
axarr[1].plot(
times,
hec_ras_record["channel_h_up"],
label="Upstream\n(HEC-RAS)",
linestyle="--",
color="darkgreen",
)
axarr[1].plot(
times,
rtc_tools_record["channel_h_dn"],
label="Downstream",
color="xkcd:dark sky blue",
)
# Format bottom axis label
axarr[-1].xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
# Shrink margins
fig.tight_layout()
# Shrink each axis by 20% and put a legend to the right of the axis
for i in range(n_subplots):
box = axarr[i].get_position()
axarr[i].set_position([box.x0, box.y0, box.width * 0.65, box.height])
axarr[i].legend(
loc="center left", bbox_to_anchor=(1, 0.5), frameon=False, prop={"size": 8}
)
plt.autoscale(enable=True, axis="x", tight=True)
# Output Plot
plt.show()
time,Inflow_Q,Level_H
2013-01-01 01:00:00,100.000000,0.000000
2013-01-01 01:15:00,100.000000,0.000000
2013-01-01 01:30:00,100.000000,0.000000
2013-01-01 01:45:00,100.000000,0.000000
2013-01-01 02:00:00,100.000000,0.000000
2013-01-01 02:15:00,100.000000,0.000000
2013-01-01 02:30:00,100.000000,0.000000
2013-01-01 02:45:00,100.000000,0.000000
2013-01-01 03:00:00,100.000000,0.000000
2013-01-01 03:15:00,100.000000,0.000000
2013-01-01 03:30:00,100.000000,0.000000
2013-01-01 03:45:00,100.000000,0.000000
2013-01-01 04:00:00,100.000000,0.000000
2013-01-01 04:15:00,100.000000,0.000000
2013-01-01 04:30:00,100.000000,0.000000
2013-01-01 04:45:00,100.000000,0.000000
2013-01-01 05:00:00,100.000000,0.000000
2013-01-01 05:15:00,100.000000,0.000000
2013-01-01 05:30:00,100.000000,0.000000
2013-01-01 05:45:00,100.000000,0.000000
2013-01-01 06:00:00,100.000000,0.000000
2013-01-01 06:15:00,100.000000,0.000000
2013-01-01 06:30:00,100.000000,0.000000
2013-01-01 06:45:00,100.000000,0.000000
2013-01-01 07:00:00,100.000000,0.000000
2013-01-01 07:15:00,100.000000,0.000000
2013-01-01 07:30:00,100.000000,0.000000
2013-01-01 07:45:00,100.000000,0.000000
2013-01-01 08:00:00,100.000000,0.000000
2013-01-01 08:15:00,100.000000,0.000000
2013-01-01 08:30:00,100.000000,0.000000
2013-01-01 08:45:00,100.000000,0.000000
2013-01-01 09:00:00,100.000000,0.000000
2013-01-01 09:15:00,100.000000,0.000000
2013-01-01 09:30:00,100.000000,0.000000
2013-01-01 09:45:00,100.000000,0.000000
2013-01-01 10:00:00,100.000000,0.000000
2013-01-01 10:15:00,100.000000,0.000000
2013-01-01 10:30:00,100.000000,0.000000
2013-01-01 10:45:00,100.000000,0.000000
2013-01-01 11:00:00,100.000000,0.000000
2013-01-01 11:15:00,100.000000,0.000000
2013-01-01 11:30:00,100.000000,0.000000
2013-01-01 11:45:00,100.000000,0.000000
2013-01-01 12:00:00,100.000000,0.000000
2013-01-01 12:15:00,100.000000,0.000000
2013-01-01 12:30:00,100.000000,0.000000
2013-01-01 12:45:00,100.000000,0.000000
2013-01-01 13:00:00,100.000000,0.000000
2013-01-01 13:15:00,100.000000,0.000000
2013-01-01 13:30:00,100.000000,0.000000
2013-01-01 13:45:00,100.000000,0.000000
2013-01-01 14:00:00,100.000000,0.000000
2013-01-01 14:15:00,100.000000,0.000000
2013-01-01 14:30:00,100.000000,0.000000
2013-01-01 14:45:00,100.000000,0.000000
2013-01-01 15:00:00,100.000000,0.000000
2013-01-01 15:15:00,100.000000,0.000000
2013-01-01 15:30:00,100.000000,0.000000
2013-01-01 15:45:00,100.000000,0.000000
2013-01-01 16:00:00,100.000000,0.000000
2013-01-01 16:15:00,100.000000,0.000000
2013-01-01 16:30:00,100.000000,0.000000
2013-01-01 16:45:00,100.000000,0.000000
2013-01-01 17:00:00,100.000000,0.000000
2013-01-01 17:15:00,100.000000,0.000000
2013-01-01 17:30:00,100.000000,0.000000
2013-01-01 17:45:00,100.000000,0.000000
2013-01-01 18:00:00,100.000000,0.000000
2013-01-01 18:15:00,100.000000,0.000000
2013-01-01 18:30:00,100.000000,0.000000
2013-01-01 18:45:00,100.000000,0.000000
2013-01-01 19:00:00,100.000000,0.000000
2013-01-01 19:15:00,100.000000,0.000000
2013-01-01 19:30:00,100.000000,0.000000
2013-01-01 19:45:00,100.000000,0.000000
2013-01-01 20:00:00,100.000000,0.000000
2013-01-01 20:15:00,100.000000,0.000000
2013-01-01 20:30:00,100.000000,0.000000
2013-01-01 20:45:00,100.000000,0.000000
2013-01-01 21:00:00,500.000000,0.000000
2013-01-01 21:15:00,500.000000,0.000000
2013-01-01 21:30:00,500.000000,0.000000
2013-01-01 21:45:00,500.000000,0.000000
2013-01-01 22:00:00,500.000000,0.000000
2013-01-01 22:15:00,500.000000,0.000000
2013-01-01 22:30:00,500.000000,0.000000
2013-01-01 22:45:00,500.000000,0.000000
2013-01-01 23:00:00,500.000000,0.000000
2013-01-01 23:15:00,500.000000,0.000000
2013-01-01 23:30:00,500.000000,0.000000
2013-01-01 23:45:00,500.000000,0.000000
2013-01-02 00:00:00,500.000000,0.000000
2013-01-02 00:15:00,500.000000,0.000000
2013-01-02 00:30:00,500.000000,0.000000
2013-01-02 00:45:00,500.000000,0.000000
2013-01-02 01:00:00,500.000000,0.000000
2013-01-02 01:15:00,500.000000,0.000000
2013-01-02 01:30:00,500.000000,0.000000
2013-01-02 01:45:00,500.000000,0.000000
2013-01-02 02:00:00,500.000000,0.000000
2013-01-02 02:15:00,500.000000,0.000000
2013-01-02 02:30:00,500.000000,0.000000
2013-01-02 02:45:00,500.000000,0.000000
2013-01-02 03:00:00,500.000000,0.000000
2013-01-02 03:15:00,500.000000,0.000000
2013-01-02 03:30:00,500.000000,0.000000
2013-01-02 03:45:00,500.000000,0.000000
2013-01-02 04:00:00,500.000000,0.000000
2013-01-02 04:15:00,500.000000,0.000000
2013-01-02 04:30:00,500.000000,0.000000
2013-01-02 04:45:00,500.000000,0.000000
2013-01-02 05:00:00,500.000000,0.000000
2013-01-02 05:15:00,500.000000,0.000000
2013-01-02 05:30:00,500.000000,0.000000
2013-01-02 05:45:00,500.000000,0.000000
2013-01-02 06:00:00,500.000000,0.000000
2013-01-02 06:15:00,500.000000,0.000000
2013-01-02 06:30:00,500.000000,0.000000
2013-01-02 06:45:00,500.000000,0.000000
2013-01-02 07:00:00,100.000000,0.000000
2013-01-02 07:15:00,100.000000,0.000000
2013-01-02 07:30:00,100.000000,0.000000
2013-01-02 07:45:00,100.000000,0.000000
2013-01-02 08:00:00,100.000000,0.000000
2013-01-02 08:15:00,100.000000,0.000000
2013-01-02 08:30:00,100.000000,0.000000
2013-01-02 08:45:00,100.000000,0.000000
2013-01-02 09:00:00,100.000000,0.000000
2013-01-02 09:15:00,100.000000,0.000000
2013-01-02 09:30:00,100.000000,0.000000
2013-01-02 09:45:00,100.000000,0.000000
2013-01-02 10:00:00,100.000000,0.000000
2013-01-02 10:15:00,100.000000,0.000000
2013-01-02 10:30:00,100.000000,0.000000
2013-01-02 10:45:00,100.000000,0.000000
2013-01-02 11:00:00,100.000000,0.000000
2013-01-02 11:15:00,100.000000,0.000000
2013-01-02 11:30:00,100.000000,0.000000
2013-01-02 11:45:00,100.000000,0.000000
2013-01-02 12:00:00,100.000000,0.000000
2013-01-02 12:15:00,100.000000,0.000000
2013-01-02 12:30:00,100.000000,0.000000
2013-01-02 12:45:00,100.000000,0.000000
2013-01-02 13:00:00,100.000000,0.000000
2013-01-02 13:15:00,100.000000,0.000000
2013-01-02 13:30:00,100.000000,0.000000
2013-01-02 13:45:00,100.000000,0.000000
2013-01-02 14:00:00,100.000000,0.000000
2013-01-02 14:15:00,100.000000,0.000000
2013-01-02 14:30:00,100.000000,0.000000
2013-01-02 14:45:00,100.000000,0.000000
2013-01-02 15:00:00,100.000000,0.000000
2013-01-02 15:15:00,100.000000,0.000000
2013-01-02 15:30:00,100.000000,0.000000
2013-01-02 15:45:00,100.000000,0.000000
2013-01-02 16:00:00,100.000000,0.000000
2013-01-02 16:15:00,100.000000,0.000000
2013-01-02 16:30:00,100.000000,0.000000
2013-01-02 16:45:00,100.000000,0.000000
2013-01-02 17:00:00,100.000000,0.000000
2013-01-02 17:15:00,100.000000,0.000000
2013-01-02 17:30:00,100.000000,0.000000
2013-01-02 17:45:00,100.000000,0.000000
2013-01-02 18:00:00,100.000000,0.000000
2013-01-02 18:15:00,100.000000,0.000000
2013-01-02 18:30:00,100.000000,0.000000
2013-01-02 18:45:00,100.000000,0.000000
2013-01-02 19:00:00,100.000000,0.000000
2013-01-02 19:15:00,100.000000,0.000000
2013-01-02 19:30:00,100.000000,0.000000
2013-01-02 19:45:00,100.000000,0.000000
2013-01-02 20:00:00,100.000000,0.000000
2013-01-02 20:15:00,100.000000,0.000000
2013-01-02 20:30:00,100.000000,0.000000
2013-01-02 20:45:00,100.000000,0.000000
2013-01-02 21:00:00,100.000000,0.000000
2013-01-02 21:15:00,100.000000,0.000000
2013-01-02 21:30:00,100.000000,0.000000
2013-01-02 21:45:00,100.000000,0.000000
2013-01-02 22:00:00,100.000000,0.000000
2013-01-02 22:15:00,100.000000,0.000000
2013-01-02 22:30:00,100.000000,0.000000
2013-01-02 22:45:00,100.000000,0.000000
2013-01-02 23:00:00,100.000000,0.000000
2013-01-02 23:15:00,100.000000,0.000000
2013-01-02 23:30:00,100.000000,0.000000
2013-01-02 23:45:00,100.000000,0.000000
2013-01-03 00:00:00,100.000000,0.000000
2013-01-03 00:15:00,100.000000,0.000000
2013-01-03 00:30:00,100.000000,0.000000
model Example
// Elements
Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal Channel(
H_b_down = -5.0,
H_b_up = -5.0,
friction_coefficient = 0.045,
use_manning = true,
length = 10000,
theta = theta,
use_inertia = true,
use_convective_acceleration = false,
use_upwind = false,
n_level_nodes = 11,
uniform_nominal_depth = 5.0,
bottom_width_down = 30,
bottom_width_up = 30,
left_slope_angle_up = 45,
left_slope_angle_down = 45,
right_slope_angle_up = 45,
right_slope_angle_down = 45,
semi_implicit_step_size = 0.0 // Here, this does matter: step_size
) annotation(Placement(visible = true, transformation(origin = {0, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level Level annotation(Placement(visible = true, transformation(origin = {60, 8}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge Discharge annotation(Placement(visible = true, transformation(origin = {-60, 8}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
// Inputs
input Real Inflow_Q(fixed=true) = Discharge.Q;
input Real Level_H(fixed=true) = Level.H;
parameter Real theta;
parameter Real step_size;
// Output Channel states
output Real Channel_Q_up = Discharge.Q;
output Real Channel_Q_dn = Level.HQ.Q;
output Real Channel_H_up = Discharge.HQ.H;
output Real Channel_H_dn = Level.H;
equation
connect(Channel.HQDown, Level.HQ) annotation(Line(points = {{8, 0}, {60, 0}, {60, 0}, {60, 0}}, color = {0, 0, 255}));
connect(Discharge.HQ, Channel.HQUp) annotation(Line(points = {{-60, 0}, {-8, 0}, {-8, 0}, {-8, 0}}, color = {0, 0, 255}));
initial equation
Channel.Q = fill(Inflow_Q, Channel.n_level_nodes + 1);
der(Channel.H) = fill(0.0, Channel.n_level_nodes);
end Example;
This diff is collapsed.
This diff is collapsed.
from rtctools.optimization.collocated_integrated_optimization_problem import (
CollocatedIntegratedOptimizationProblem,
)
from rtctools.optimization.csv_mixin import CSVMixin
from rtctools.optimization.homotopy_mixin import HomotopyMixin
from rtctools.optimization.modelica_mixin import ModelicaMixin
from rtctools.util import run_optimization_problem
class Example(
HomotopyMixin, CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem
):
def parameters(self, ensemble_member):
p = super().parameters(ensemble_member)
times = self.times()
p['step_size'] = times[1] - times[0]
p['Channel.use_convective_acceleration'] = self.use_convective_acceleration
p['Channel.use_upwind'] = self.use_upwind
return p
def constraints(self, ensemble_member):
constraints = super().constraints(ensemble_member)
times = self.times()
# Extract the number of nodes in the channel
parameters = self.parameters(ensemble_member)
n_level_nodes = int(parameters["Channel.n_level_nodes"])
# To Mimic HEC-RAS behaviour, enforce steady state both at t0 and at t1.
for i in range(n_level_nodes):
state = f"Channel.H[{i + 1}]"
constraints.append(
(self.state_at(state, times[0]) - self.state_at(state, times[1]), 0, 0)
)
return constraints
class ExampleInertialWave(Example):
"""Inertial wave equation (no convective acceleration)"""
model_name = 'Example'
use_convective_acceleration = False
use_upwind = False
timeseries_export_basename = "timeseries_export_inertial_wave"
class ExampleSaintVenant(Example):
"""Saint Venant equation. Convective acceleration discretized with central differences"""
model_name = 'Example'
use_convective_acceleration = True
use_upwind = False
timeseries_export_basename = "timeseries_export_saint_venant"
class ExampleSaintVenantUpwind(Example):
"""Saint Venant equation. Convective acceleration discretized with upwind scheme"""
model_name = 'Example'
use_convective_acceleration = True
use_upwind = True
timeseries_export_basename = "timeseries_export_saint_venant_upwind"
run_optimization_problem(ExampleInertialWave)
run_optimization_problem(ExampleSaintVenant)
run_optimization_problem(ExampleSaintVenantUpwind)
......@@ -44,7 +44,8 @@ setup(
install_requires = ["casadi >= 3.2.0",
"numpy >= 1.14.0",
"scipy >= 1.0.0",
"pymoca >= 0.4.1"],
"pymoca >= 0.4.1",
"rtc-tools-channel-flow >= 1.1.0a1"],
tests_require = ['pytest', 'pytest-runner'],
python_requires='>=3.5',
cmdclass = versioneer.get_cmdclass(),
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment