Merewether Flood Case Study Example

Here we look at a case study of a flood in the community of Merewether near Newcastle NSW. We will add a flow using an Inlet_operator and extract flow details at various points by interagating the sww file which is produced during the ANUGA run.

This example is based on the the validation test merewether, provided in the ANUGA distribution.

The data for this example is available from our anuga-clinic repository.

Setup Notebook for Visualisation and Animation

We are using the format of a jupyter notebook. As such we need to setup inline matplotlib plotting and animation.

import numpy as np
import os
import matplotlib.pyplot as plt

%matplotlib inline

# Allow inline jshtml animations
from matplotlib import rc
rc('animation', html='jshtml')

Import ANUGA

We assume that anuga has been installed. If so we can import anuga.

import anuga
import pathlib

# Root of the anuga_core repository (works for editable installs and Codespaces)
ANUGA_CORE = str(pathlib.Path(anuga.__file__).parent.parent)

data_dir = os.path.join(ANUGA_CORE, 'examples', 'data', 'merewether')
print(f"Data directory: {data_dir}")

Data directory: /home/steve/anuga_core/examples/data/merewether

Read in Data

We have copied over some topography data and extent data from our anuga-clinic notebook repository.

Let’s read that in and create a mesh associated with it.

# Polygon defining broad area of interest
bounding_polygon = anuga.read_polygon(anuga.join(data_dir,'extent.csv'))


# Polygon defining particular area of interest
merewether_polygon = anuga.read_polygon(anuga.join(data_dir,'merewether.csv'))


# Elevation Data
topography_file = anuga.join(data_dir,'topography1.asc')


# Resolution for most of the mesh
base_resolution = 80.0  # m^2

# Resolution in particular area of interest
merewether_resolution = 25.0 # m^2

interior_regions = [[merewether_polygon, merewether_resolution]]

Create and View Domain

Note that we use a base_resolution to ensure a reasonable refinement over the whole region, and we use interior_regions to refine the mesh in the area of interest. In this case we pass a list of polygon, resolution pairs.

domain = anuga.create_domain_from_regions(
            bounding_polygon,
            boundary_tags={'south': [0],
                           'east':  [1],
                           'north':    [2],
                           'west':   [3]},
            maximum_triangle_area=base_resolution,
            interior_regions=interior_regions)


domain.set_name('merewether1') # Name of sww file
domain.set_epsg(32756) # geo_reference, could set UTM zone and hemisphere

domain.triplot(linewidth = 0.4)

(<Figure size 640x480 with 1 Axes>,
 <Axes: xlabel='Easting (m)', ylabel='Northing (m)'>,
 [<matplotlib.lines.Line2D at 0x719e4f1c5a90>,
  <matplotlib.lines.Line2D at 0x719e4f1c5be0>])

Setup Initial Conditions

We have to setup the values of various quantities associated with the domain. In particular we need to setup the elevation the elevation of the bed or the bathymetry. In this case we will do this using the DEM file topography1.asc .

domain.set_quantity('elevation', filename=topography_file, location='centroids') # Use function for elevation
domain.set_quantity('friction', 0.01, location='centroids')                        # Constant friction
domain.set_quantity('stage', expression='elevation', location='centroids')         # Dry Bed

fig, ax, im = domain.tripcolor(facecolors=domain.elev, cmap='Greys_r')
plt.colorbar(im, ax=ax)
plt.title("Elevation")

Text(0.5, 1.0, 'Elevation')

Setup Boundary Conditions

The rectangular domain has 4 tagged boundaries, left, top, right and bottom. We need to set boundary conditons for each of these tagged boundaries. We can set Transmissive type BC on the outflow boundaries and reflective on the others.

Br = anuga.Reflective_boundary(domain)
Bt = anuga.Transmissive_boundary(domain)

domain.set_boundary({'south':   Br,
                     'east':    Bt, # outflow
                     'north':   Bt, # outflow
                     'west':    Br})

Setup Inflow

We need some water to flow. The easiest way to input a specified amount of water is via an Inlet_operator where we can specify a discharge Q.

# Setup inlet flow
center = (382270.0, 6354285.0)
radius = 10.0
region0 = anuga.Region(domain, center=center, radius=radius)
fixed_inflow = anuga.Inlet_operator(domain, region0 , Q=19.7)

Run the Evolution

We evolve using a for statement, which evolves the quantities using the shallow water wave solver. The calculation yields every yieldstep seconds, up to a given duration.

for t in domain.evolve(yieldstep=20, duration=300):

    #dplotter.plot_depth_frame()
    domain.save_depth_frame(vmin=0.0, vmax=1.0)

    domain.print_timestepping_statistics()


# Read in the png files stored during the evolve loop
domain.make_depth_animation()

Time = 320.0000 (sec), delta t in [0.14940474, 0.14942357] (s), steps=134 (0s), elapsed (0s), eta (0s), mem=259MB
Time = 340.0000 (sec), delta t in [0.14938206, 0.14940414] (s), steps=134 (0s), elapsed (0s), eta (0s), mem=261MB
Time = 360.0000 (sec), delta t in [0.14937154, 0.14938394] (s), steps=134 (0s), elapsed (0s), eta (0s), mem=265MB
Time = 380.0000 (sec), delta t in [0.14936630, 0.14937438] (s), steps=134 (0s), elapsed (1s), eta (0s), mem=267MB
Time = 400.0000 (sec), delta t in [0.14936324, 0.14937220] (s), steps=134 (0s), elapsed (1s), eta (0s), mem=270MB
Time = 420.0000 (sec), delta t in [0.14936031, 0.14936849] (s), steps=134 (0s), elapsed (1s), eta (0s), mem=273MB
Time = 440.0000 (sec), delta t in [0.14935862, 0.14936705] (s), steps=134 (0s), elapsed (2s), eta (0s), mem=276MB
Time = 460.0000 (sec), delta t in [0.14933473, 0.14937451] (s), steps=134 (0s), elapsed (2s), eta (0s), mem=278MB
Time = 480.0000 (sec), delta t in [0.14936045, 0.14938262] (s), steps=134 (0s), elapsed (2s), eta (0s), mem=281MB
Time = 500.0000 (sec), delta t in [0.14935164, 0.14938343] (s), steps=134 (0s), elapsed (3s), eta (0s), mem=283MB
Time = 520.0000 (sec), delta t in [0.14933428, 0.14935610] (s), steps=134 (0s), elapsed (3s), eta (0s), mem=285MB
Time = 540.0000 (sec), delta t in [0.14932554, 0.14934195] (s), steps=134 (0s), elapsed (4s), eta (0s), mem=287MB
Time = 560.0000 (sec), delta t in [0.14932088, 0.14933516] (s), steps=134 (0s), elapsed (4s), eta (0s), mem=288MB
Time = 580.0000 (sec), delta t in [0.14931833, 0.14933023] (s), steps=134 (0s), elapsed (4s), eta (0s), mem=290MB
Time = 600.0000 (sec), delta t in [0.14931842, 0.14932949] (s), steps=134 (0s), elapsed (5s), eta (0s), mem=291MB

Once Loop Reflect

SWW File

The evolve loop saves the quantites at the end of each yield step to an sww file, with name domain name + extension sww. In this case the sww file is merewether1.sww.

An sww file can be viewed via our 3D anuga-viewer application, via the crayfish plugin for QGIS, or simply read back into python using netcdf commands.

For this clinic we have provided a wrapper called an SWW_plotter to provide easy acces to the saved quantities, stage, elev, depth, xmom, ymom, xvel, yvel, speed which are all time slices of centroid values, and a time variable.

# Create a wrapper for contents of sww file
swwfile = 'merewether1.sww'
splotter = anuga.SWW_plotter(swwfile)


# Plot Depth and Speed at the last time slice
plt.subplot(121)
splotter.triang.set_mask(None)
plt.tripcolor(splotter.triang,
              facecolors = splotter.depth[-1,:],
              cmap='viridis')

plt.title("Depth")


plt.subplot(122)
splotter.triang.set_mask(None)
plt.tripcolor(splotter.triang,
              facecolors = splotter.speed[-1,:],
              cmap='viridis')

plt.title("Speed");

Comparison

The data file ObservationPoints.csv contains some comparison depth data from Australian Rain and Runoff. Let’s plot the depth for our simulation against the comparison data.

point_observations = np.genfromtxt(
    os.path.join(data_dir,'ObservationPoints.csv'),
    delimiter=",",skip_header=1)

# Convert to absolute corrdinates
xc = splotter.xc + splotter.xllcorner
yc = splotter.yc + splotter.yllcorner

nearest_points = []
for row in point_observations:
    nearest_points.append(np.argmin( (xc-row[0])**2 + (yc-row[1])**2 ))

loc_id = point_observations[:,2]

fig, ax = plt.subplots()
ax.plot(loc_id, point_observations[:,4], '*r', label='ARR')
ax.plot(loc_id, point_observations[:,5], '*b', label='Tuflow')
ax.plot(loc_id, splotter.stage[-1,nearest_points], '*g', label='Anuga')

plt.xticks(range(0,5))
plt.xlabel('ID')
plt.ylabel('Stage')
ax.legend()

plt.show()

Flow with Houses

We have polygonal data which specifies the location of a number of structures (homes) in our study. We can consider the flow in which those houses are cut out of the simulation.

First we read in the house polygonal data. To maintain a small mesh size we will only read in structures with an area grester than 60 m^2.

# Read in house polygons from data directory and retain those of area > 60 m^2

import glob
house_files = glob.glob(os.path.join(data_dir,'house*.csv'))

house_polygons = []
for hf in house_files:
  house_poly = anuga.read_polygon(hf)
  poly_area = anuga.polygon_area(house_poly)

  # Leave out some small houses
  if poly_area > 60:
    house_polygons.append(house_poly)

Create Domain

To incorporate the housing information, we will cutout the polygons representing the houses. This is done by passing the list of house polygons to the interior_holes argument of the anuga.create_domain_from_regions procedure.

This will produce a new tagged boundary region called interior. We will have to assign a boundsry condition to this new boundary region.

# Resolution for most of the mesh
base_resolution = 20.0  # m^2

# Resolution in particular area of interest
merewether_resolution = 10.0 # m^2

domain = anuga.create_domain_from_regions(
            bounding_polygon,
            boundary_tags={'bottom': [0],
                           'right':  [1],
                           'top':    [2],
                           'left':   [3]},
            maximum_triangle_area=base_resolution,
            interior_holes=house_polygons,
            interior_regions=interior_regions)


domain.set_name('merewether2') # Name of sww file
domain.set_epsg(32756) # geo_reference, could set UTM zone and hemisphere
domain.set_low_froude(1)

# Setup Initial Conditions
domain.set_quantity('elevation', filename=topography_file, location='centroids') # Use function for elevation
domain.set_quantity('friction', 0.01, location='centroids')                        # Constant friction
domain.set_quantity('stage', expression='elevation', location='centroids')         # Dry Bed

# Setup BC
Br = anuga.Reflective_boundary(domain)
Bt = anuga.Transmissive_boundary(domain)


# NOTE: We need to assign a BC to the interior boundary region.
domain.set_boundary({'bottom':   Br,
                     'right':    Bt, # outflow
                     'top':      Bt, # outflow
                     'left':     Br,
                     'interior': Br})


# Setup inlet flow
center = (382270.0, 6354285.0)
radius = 10.0
region0 = anuga.Region(domain, center=center, radius=radius)
fixed_inflow = anuga.Inlet_operator(domain, region0 , Q=19.7)


domain.triplot(linewidth = 0.4)

(<Figure size 640x480 with 1 Axes>,
 <Axes: xlabel='Easting (m)', ylabel='Northing (m)'>,
 [<matplotlib.lines.Line2D at 0x719e3ff98050>,
  <matplotlib.lines.Line2D at 0x719e3ff981a0>])

Evolve

for t in domain.evolve(yieldstep=20, duration=300):

    #dplotter.plot_depth_frame()
    domain.save_depth_frame(vmin=0.0, vmax=1.0)

    domain.print_timestepping_statistics()


# Read in the png files stored during the evolve loop
domain.make_depth_animation()

Time = 320.0000 (sec), delta t in [0.09267182, 0.09269317] (s), steps=216 (0s), elapsed (0s), eta (0s), mem=344MB
Time = 340.0000 (sec), delta t in [0.09266801, 0.09270692] (s), steps=216 (0s), elapsed (1s), eta (1s), mem=344MB
Time = 360.0000 (sec), delta t in [0.09266763, 0.09270886] (s), steps=216 (0s), elapsed (2s), eta (1s), mem=344MB
Time = 380.0000 (sec), delta t in [0.09266775, 0.09270829] (s), steps=216 (0s), elapsed (2s), eta (1s), mem=344MB
Time = 400.0000 (sec), delta t in [0.09266779, 0.09270732] (s), steps=216 (0s), elapsed (3s), eta (1s), mem=344MB
Time = 420.0000 (sec), delta t in [0.09266794, 0.09270683] (s), steps=216 (0s), elapsed (4s), eta (1s), mem=344MB
Time = 440.0000 (sec), delta t in [0.09266793, 0.09270668] (s), steps=216 (1s), elapsed (5s), eta (1s), mem=344MB
Time = 460.0000 (sec), delta t in [0.09266809, 0.09270657] (s), steps=216 (0s), elapsed (5s), eta (1s), mem=344MB
Time = 480.0000 (sec), delta t in [0.09266806, 0.09270653] (s), steps=216 (0s), elapsed (6s), eta (1s), mem=344MB
Time = 500.0000 (sec), delta t in [0.09266818, 0.09270654] (s), steps=216 (0s), elapsed (7s), eta (1s), mem=344MB
Time = 520.0000 (sec), delta t in [0.09266814, 0.09270654] (s), steps=216 (0s), elapsed (7s), eta (1s), mem=344MB
Time = 540.0000 (sec), delta t in [0.09266824, 0.09270658] (s), steps=216 (0s), elapsed (8s), eta (0s), mem=344MB
Time = 560.0000 (sec), delta t in [0.09266820, 0.09270656] (s), steps=216 (0s), elapsed (9s), eta (0s), mem=344MB
Time = 580.0000 (sec), delta t in [0.09266829, 0.09270663] (s), steps=216 (0s), elapsed (9s), eta (0s), mem=344MB
Time = 600.0000 (sec), delta t in [0.09266824, 0.09270658] (s), steps=216 (0s), elapsed (10s), eta (0s), mem=344MB

Once Loop Reflect

Read in SWW File and Compare

Perhaps not conclusive, but with the houses the anuga results, especially for id point 0, are much closer to the comparison results. Note that we are running with a very coarse mesh for this case study.

# Create a wrapper for contents of sww file
swwfile2 = 'merewether2.sww'
splotter2 = anuga.SWW_plotter(swwfile2)

# Convert to absolute corrdinates
xc = splotter2.xc + splotter2.xllcorner
yc = splotter2.yc + splotter2.yllcorner

nearest_points_2 = []
for row in point_observations:
    nearest_points_2.append(np.argmin( (xc-row[0])**2 + (yc-row[1])**2 ))

loc_id = point_observations[:,2]

fig, ax = plt.subplots()
ax.plot(loc_id, point_observations[:,4], '*r', label='ARR')
ax.plot(loc_id, point_observations[:,5], '*b', label='Tuflow')
ax.plot(loc_id, splotter2.stage[-1,nearest_points_2], '*g', label='Anuga1')
ax.plot(loc_id, splotter.stage[-1,nearest_points], '*k', label='Anuga0')


plt.xticks(range(0,5))
plt.xlabel('ID')
plt.ylabel('Stage')
ax.legend()

plt.show()