"""
Author: HECE - University of Liege, Stéphane Champailler, Pierre Archambeau
Date: 2024
Copyright (c) 2024 University of Liege. All rights reserved.
This script and its content are protected by copyright law. Unauthorized
copying or distribution of this file, via any medium, is strictly prohibited.
"""
from enum import Enum
import numpy as np
from .simple_simulation import SimulationDuration, TimeStepStrategy, SimpleSimulation, BoundaryConditionsTypes
[docs]
class WallDirection(Enum):
""" Walls follow x-axis """
""" Walls follow y-axis """
[docs]
HORIZONTAL_AND_VERTICAL = 3
""" Walls follow x and y-axis """
[docs]
def make_closed_pool(width, height, pos_x, pos_y, pool_width, pool_height, wall_height=10,
mode = "cube") -> SimpleSimulation:
""" Make a cube of water, surrounded by empty space, surrounded by walls,
surrounded by empty space.
The problem will be pool_width by pool_height but it will be placed
onto a grid (width*height) at (pos_x,pos_y)
"""
assert width >= 6*2 and height >= 6*2, "sim must be wide enough to fit everything inside it"
assert pool_width <= width and pool_height <= height, "problem size doesn't fit inside the grid"
# Make the cube of water surrounded by a wall.
shape = (pool_width, pool_height)
# Surround the swimming pool by a wall, 2 cells wide
bathymetry = np.zeros(shape, np.float32)
bathymetry[ 2:-2, 2:-2 ] = wall_height
bathymetry[ 4:-4, 4:-4 ] = 0
# And now, we drop the small pool inside the big simulation
sim = SimpleSimulation(width, height)
sim.param_dx = 1.0
sim.param_dy = 1.0
sim.h = np.zeros(sim.shape, np.float32)
sim.bathymetry = np.zeros(sim.shape, np.float32)
sim.nap = np.zeros(sim.shape, np.uint8)
sim.qx = np.zeros(sim.shape, np.float32)
sim.qy = np.zeros(sim.shape, np.float32)
if mode == "cube":
h = np.zeros(shape, np.float32)
cube_width = pool_width // 4
cube_height = pool_height // 3
assert cube_width >= 1, "widths are too small, can't build the cube of water."
assert cube_height >= 1, "heights are too small, can't build the cube of water."
h[pool_width//2 -cube_width//2 :pool_width//2 +cube_width//2,
pool_height//2-cube_height//2:pool_height//2+cube_height//2] = 1.0
elif mode == "filled":
h = np.zeros(shape, np.float32)
h[4:-4, 4:-4] = 1.0
sim.h[pos_x:pos_x+pool_width, pos_y:pos_y+pool_height] = h
sim.bathymetry[pos_x:pos_x+pool_width, pos_y:pos_y+pool_height] = bathymetry
sim.nap[pos_x:pos_x+pool_width, pos_y:pos_y+pool_height] = 1
sim.manning = np.full(sim.shape, 0.04, np.float32)
sim.nap[:2,:] = 0
sim.nap[-2:,:] = 0
sim.nap[:,:2] = 0
sim.nap[:,-2:] = 0
return sim
[docs]
def make_cube_drop(width, height, water_cube_height=10.0) -> SimpleSimulation:
""" Create a simulation with a cube of water in the middle of the domain.
- The block's size is 1/4 of the domain's size.
- The spatial resolution is 1.0 meter per cell.
- We will compute 10_000 time steps.
- No discharges at initial time.
- Topography is flat (0. m), except for the boundaries where we have a wall (100 m).
- Manning coefficient is 0.4 everywhere (very rough).
:attention: USED IN BENCHMARKING
:param width, height: dimensions of the simulation domain.
:param water_cube_height: height of the cube of water (meters).
"""
sim = SimpleSimulation(width, height)
sim.param_dx, sim.param_dy = 1.0, 1.0
sim._param_nx, sim._param_ny = width, height
sim.param_duration = SimulationDuration.from_steps(10000)
shape = (sim._param_nx, sim._param_ny)
sim.nap = np.zeros(shape, dtype=np.uint8)
sim.nap[ 2:-2, 2:-2] = True
sim._h = np.ones(shape, dtype=np.float32) * 0.0
mx = sim._param_nx // 2
nx = sim._param_nx // 4
my = sim._param_ny // 2
ny = sim._param_ny // 4
sim._h[ mx-nx:mx+nx, my-ny:my+ny] = water_cube_height
sim._qx = np.zeros(shape, dtype=np.float32)
sim._qy = np.zeros(shape, dtype=np.float32)
sim._bathymetry = np.zeros(shape, dtype=np.float32)
sim._bathymetry[:, 0:2] = 100
sim._bathymetry[:, sim._param_ny-2:sim._param_ny] = 100
sim._bathymetry[0:2, :] = 100
sim._bathymetry[sim._param_nx-2:sim._param_nx, :] = 100
sim.manning = np.ones(shape, dtype=np.float32)*0.4
return sim
[docs]
def make_swimming_pool(width:int, height:int,
water_depth:float,
wall: WallDirection,
manning:float= 0.4,
dx:float= 1, dy:float= 1) -> SimpleSimulation:
""" Create a rectangular swimming pool surrounded by walls.
:Note:
- the meshes at x,y=0 and x,y=W/H-1 are not used.
- the meshes at x,y=1 and x,y=W/H-2 are filled with walls.
Therefore the smallest fully enclosed swimming pool one can build is 5x5 meshes.
Manning coefficient is set to 0.4 everywhere.
:param width: dimensions of the swimming pool along X.
:param height: dimensions of the swimming pool along Y.
:param water_depth: the swimming pool will be filled with that height of water (meters).
:param wall: which walls should be build around the pool ?
You can choose to build only those along x-axis (hoorizontal),
y-axis (vertical) or both.
:param manning: the Manning coefficient of the water in the pool (default = 0.4).
"""
assert not (wall is None and water_depth > 0), "If no walls, then no water"
sim = SimpleSimulation(width, height)
sim._param_nx, sim._param_ny = width, height
sim.param_dx, sim.param_dy = float(dx), float(dy)
shape = (sim._param_nx, sim._param_ny)
sim.nap = np.zeros(shape, dtype=np.uint8)
sim.nap[ 1:sim._param_nx-1, 1:sim._param_ny-1] = 1
sim._h = np.zeros(shape, dtype=np.float32)
sim._qx = np.zeros(shape, dtype=np.float32)
sim._qy = np.zeros(shape, dtype=np.float32)
sim._bathymetry = np.zeros(shape, dtype=np.float32)
if wall == WallDirection.HORIZONTAL:
sim._h[1:sim._param_nx-1, 2:sim._param_ny-2] = water_depth
elif wall == WallDirection.VERTICAL:
sim._h[2:sim._param_nx-2, 1:sim._param_ny-1] = water_depth
elif wall == WallDirection.HORIZONTAL_AND_VERTICAL:
sim._h[2:sim._param_nx-2, 2:sim._param_ny-2] = water_depth
if wall in (WallDirection.HORIZONTAL, WallDirection.HORIZONTAL_AND_VERTICAL):
sim._bathymetry[:, 0:2] = 100
sim._bathymetry[:, sim._param_ny-2:sim._param_ny] = 100
if wall in (WallDirection.VERTICAL, WallDirection.HORIZONTAL_AND_VERTICAL):
sim._bathymetry[0:2, :] = 100
sim._bathymetry[sim._param_nx-2:sim._param_nx, :] = 100
sim._manning = np.ones(shape, dtype=np.float32)*manning
return sim
[docs]
def make_channel(width:int, height:int,
water_depth:float,
direction: WallDirection,
manning:float= 0.04,
dx:float= 1, dy:float= 1) -> SimpleSimulation:
""" Create a channel of water surrounded by walls.
:note:
- one cell all around the domain is not computed.
- the walls are built on the second cell from the border.
- the water is put inside the walls.
Therefore the smallest fully enclosed channel one can build is 5 cells wide.
:param width: dimensions of the channel along X.
:param height: dimensions of the channel along Y.
:param water_depth: the channel will be filled with
that height of water (meters).
:param direction: which walls should be build around the channel ?
You can choose to build only those along x-axis (hoorizontal),
y-axis (vertical) or both.
"""
assert not (direction is None and water_depth > 0), "If no walls, then no water"
assert direction in (WallDirection.HORIZONTAL, WallDirection.VERTICAL), "Only horizontal or vertical channels are supported"
sim = make_swimming_pool(width, height, water_depth, direction, manning, dx, dy)
return sim
[docs]
def add_triangular_bridge2channel(sim: SimpleSimulation, pos:float = None,
width:float= 20., roof_elevation_min:float = 5., roof_elevation_max:float = 6.):
""" Add a triangular bridge to the simulation.
:param sim: the simulation to modify.
:param pos: position of the bridge (in meters). If None, the bridge is placed at the center of the simulation.
:param width: width of the bridge (in meters).
:param elevation_min: minimum elevation of the bridge's roof (in meters) at the center.
:param elevation_max: maximum elevation of the bridge's roof (in meters) at the extrema.
"""
if sim.param_nx > sim.param_ny:
direction = WallDirection.HORIZONTAL
else:
direction = WallDirection.VERTICAL
sim.make_void_bridge_roof()
if pos is None:
if direction == WallDirection.HORIZONTAL:
pos = sim.param_nx // 2
start_i = pos - int(width // sim.param_dx) // 2
end_i = pos + int(width // sim.param_dx) // 2
start_j = 2
end_j = sim.param_ny - 2
#ensure that the bridge is an odd lenght cells wide
if (end_i - start_i) % 2 == 0:
start_i -= 1
center_i = (start_i + end_i) // 2
for j in range(start_j, end_j):
sim._bridge_roof_elevation[start_i:end_i, j] = np.concatenate([np.linspace(roof_elevation_max, roof_elevation_min, center_i+1 - start_i),
np.linspace(roof_elevation_min, roof_elevation_max, end_i - center_i)[1:]])
else:
pos = sim.param_ny // 2
start_i = 2
end_i = sim.param_nx - 2
start_j = pos - int(width // sim.param_dy) // 2
end_j = pos + int(width // sim.param_dy) // 2
#ensure that the bridge is an odd lenght cells wide
if (end_j - start_j) % 2 == 0:
start_j -= 1
center_j = (start_j + end_j) // 2
for i in range(start_i, end_i):
sim._bridge_roof_elevation[i, start_j:end_j] = np.concatenate([np.linspace(roof_elevation_max, roof_elevation_min, center_j+1 - start_j),
np.linspace(roof_elevation_min, roof_elevation_max, end_j - center_j)[1:]])
sim.convert_roof_array2params()
[docs]
def add_arch_bridge2channel(sim: SimpleSimulation, pos:float = None,
width:float= 20., roof_elevation_min:float = 5., roof_elevation_max:float = 6.):
""" Add an single arch bridge to the simulation.
:param sim: the simulation to modify.
:param pos: position of the bridge (in meters). If None, the bridge is placed at the center of the simulation.
:param width: width of the bridge (in meters).
:param elevation_min: minimum elevation of the bridge's roof (in meters) at the extrema.
:param elevation_max: maximum elevation of the bridge's roof (in meters) at the center.
"""
if sim.param_nx > sim.param_ny:
direction = WallDirection.HORIZONTAL
else:
direction = WallDirection.VERTICAL
sim.make_void_bridge_roof()
if pos is None:
if direction == WallDirection.HORIZONTAL:
pos = sim.param_nx // 2
start_i = pos - int(width // sim.param_dx) // 2
end_i = pos + int(width // sim.param_dx) // 2
# arch bridge as polynomial with roof_elevation_min at extremes and horizontal tangent at center
# z = a * (x - x0) ** 2 + b
# b = roof_elevation_max
center_j = sim.param_ny // 2
start_j = 2
end_j = sim.param_ny - 2
a = (roof_elevation_min - roof_elevation_max) / (center_j - start_j) ** 2
b = roof_elevation_max
z = np.asarray([a * (j - center_j) ** 2 + b for j in range(start_j, end_j)])
for i in range(start_i, end_i):
sim._bridge_roof_elevation[i, 2:-2] = z
else:
pos = sim.param_ny // 2
start_j = pos - int(width // sim.param_dy) // 2
end_j = pos + int(width // sim.param_dy) // 2
# arch bridge as polynomial with roof_elevation_min at extremes and horizontal tangent at center
# z = a * (x - x0) ** 2 + b
# b = roof_elevation_max
center_i = sim.param_nx // 2
start_i = 2
end_i = sim.param_nx - 2
a = (roof_elevation_min - roof_elevation_max) / (center_i - start_i) ** 2
b = roof_elevation_max
z = np.asarray([a * (i - center_i) ** 2 + b for i in range(start_i, end_i)])
for j in range(start_j, end_j):
sim._bridge_roof_elevation[2:-2, j] = z
sim.convert_roof_array2params()
[docs]
def add_multi_arch_bridge2channel(sim: SimpleSimulation, pos:float = None,
width:float= 20., nb_arches:int = 2,
roof_elevation_min:float = 5., roof_elevation_max:float = 6.):
""" Add a multi-arch bridge to the simulation.
:param sim: the simulation to modify.
:param pos: position of the bridge (in meters). If None, the bridge is placed at the center of the simulation.
:param width: width of the bridge (in meters).
:param nb_arches: number of arches in the bridge.
:param elevation_min: minimum elevation of the bridge's roof (in meters) at the extrema.
:param elevation_max: maximum elevation of the bridge's roof (in meters) at the center.
"""
if sim.param_nx > sim.param_ny:
direction = WallDirection.HORIZONTAL
else:
direction = WallDirection.VERTICAL
sim.make_void_bridge_roof()
if pos is None:
if direction == WallDirection.HORIZONTAL:
pos = sim.param_nx // 2
start_i = pos - int(width // sim.param_dx) // 2
end_i = pos + int(width // sim.param_dx) // 2
# arch bridge as polynomial with roof_elevation_min at extremes and horizontal tangent at center
# z = a * (x - x0) ** 2 + b
# b = roof_elevation_max
z = []
w = sim.param_ny - 4
center_j = w // nb_arches // 2 + 2
start_j = 2
end_j = start_j + w // nb_arches
for curarch in range(nb_arches):
a = (roof_elevation_min - roof_elevation_max) / (center_j - start_j) ** 2
b = roof_elevation_max
z.append(np.asarray([a * (j - center_j) ** 2 + b for j in range(start_j, end_j)]))
center_j += w // nb_arches
start_j += w // nb_arches
end_j += w // nb_arches
z = np.concatenate(z)
for i in range(start_i, end_i):
sim._bridge_roof_elevation[i, 2:-2] = z
else:
pos = sim.param_ny // 2
start_j = pos - int(width // sim.param_dy) // 2
end_j = pos + int(width // sim.param_dy) // 2
# arch bridge as polynomial with roof_elevation_min at extremes and horizontal tangent at center
# z = a * (x - x0) ** 2 + b
# b = roof_elevation_max
z = []
w = sim.param_nx - 4
center_i = w // nb_arches // 2 + 2
start_i = 2
end_i = start_i + w // nb_arches
for curarch in range(nb_arches):
a = (roof_elevation_min - roof_elevation_max) / (center_i - start_i) ** 2
b = roof_elevation_max
z.append(np.asarray([a * (i - center_i) ** 2 + b for i in range(start_i, end_i)]))
center_i += w // nb_arches
start_i += w // nb_arches
end_i += w // nb_arches
z = np.concatenate(z)
for j in range(start_j, end_j):
sim._bridge_roof_elevation[2:-2, j] = z
sim.convert_roof_array2params()