How to use anuga for simulating the downstream impacts of discharging water from a reservoir
chooron opened this issue · 5 comments
I am studying an area with a reservoir and want to investigate the changes in runoff when the reservoir discharges water downstream. However, if I set the inflow to a low value, the water will not flow downstream, and how to set the direction of the inflow
Here is my code:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import anuga
import os
from anuga.file.sww import SWW_file
from anuga.utilities import animate
from matplotlib import rc
rc('animation', html='jshtml')
# set the dem file path
asc_file_path = r'../../data/asc/lintao_dem_90.asc'
file = open(asc_file_path)
polygon_dict = {}
for i in range(5):
temp_str = file.readline()
temp_str_list = temp_str.split(' ')
polygon_dict[temp_str_list[0]] = float(temp_str_list[-1].split('\n')[0])
# set the bounding polygon
bounding_polygon = [(polygon_dict['xllcorner'], polygon_dict['yllcorner']),
(polygon_dict['xllcorner'] + polygon_dict['ncols'] * polygon_dict['cellsize'],
polygon_dict['yllcorner']),
(polygon_dict['xllcorner'] + polygon_dict['ncols'] * polygon_dict['cellsize'],
polygon_dict['yllcorner'] + polygon_dict['nrows'] * polygon_dict['cellsize']),
(polygon_dict['xllcorner'],
polygon_dict['yllcorner'] + polygon_dict['nrows'] * polygon_dict['cellsize'])]
# Resolution for most of the mesh
scale = 50
base_resolution = polygon_dict['cellsize'] ** 2 * scale
# build domain
domain = anuga.create_domain_from_regions(
bounding_polygon,
boundary_tags={'bottom': [0],
'right': [1],
'top': [2],
'left': [3]},
use_cache=False,
maximum_triangle_area=base_resolution,
mesh_filename='cache/lintao_mesh.msh'
# interior_regions=interior_regions
)
# visual the sww file
domain.set_name('lintao') # Name of sww file
domain.set_datadir('cache')
dplotter = animate.Domain_plotter(domain)
plt.triplot(dplotter.triang, linewidth=0.4)
plt.show()
domain.set_quantity('elevation', filename=asc_file_path, location='centroids') # Use function for elevation
domain.set_quantity('friction', 0.01, location='centroids') # Constant friction
domain.set_quantity('stage', 5.0) # Dry Bed
# anuga.Set_stage(domain, stage=500)
plt.tripcolor(dplotter.triang,
facecolors=dplotter.elev,
cmap='Greys_r')
plt.colorbar()
plt.title("Elevation")
plt.show()
Br = anuga.Reflective_boundary(domain)
Bt = anuga.Transmissive_boundary(domain)
domain.set_boundary({'bottom': Bt,
'right': Br, # outflow
'top': Br, # outflow
'left': Bt})
# Setup inlet flow
center = (polygon_dict['xllcorner'] + 37000, polygon_dict['yllcorner'] + 22000)
radius = polygon_dict['cellsize'] * 50
# region0 = anuga.Region(domain, center=center, radius=radius)
region0 = anuga.Region(domain, polygon=[
(polygon_dict['xllcorner'] + 31000, polygon_dict['yllcorner'] + 19000),
(polygon_dict['xllcorner'] + 31700, polygon_dict['yllcorner'] + 19000),
(polygon_dict['xllcorner'] + 31700, polygon_dict['yllcorner'] + 20000),
(polygon_dict['xllcorner'] + 31000, polygon_dict['yllcorner'] + 20000),
])
fixed_inflow = anuga.Inlet_operator(domain, [[polygon_dict['xllcorner'] + 31000, polygon_dict['yllcorner'] + 19000],
[polygon_dict['xllcorner'] + 31700, polygon_dict['yllcorner'] + 20000]],
Q=lambda t: 3e3)
gauge = [[polygon_dict['xllcorner'] + 35000, polygon_dict['yllcorner'] + 22000],
[polygon_dict['xllcorner'] + 35000, polygon_dict['yllcorner'] + 23000],
[polygon_dict['xllcorner'] + 35000, polygon_dict['yllcorner'] + 23000]]
g1_id = domain.get_triangle_containing_point(gauge[0])
g2_id = domain.get_triangle_containing_point(gauge[1])
g3_id = domain.get_triangle_containing_point(gauge[2])
domain.set_store(True)
for t in domain.evolve(yieldstep=50, duration=1000):
# dplotter.plot_depth_frame()
dplotter.save_depth_frame(vmin=0.0, vmax=20)
domain.print_timestepping_statistics()
# Read in the png files stored during the evolve loop
swwplotter = animate.SWW_plotter("lintao.sww")
triang = swwplotter.triang
depth = swwplotter.depth
swwplotter.triang.set_mask(None)
plt.tripcolor(swwplotter.triang,
facecolors=swwplotter.depth[-1, :],
cmap='viridis')
plt.title("Depth")
plt.show()
plt.plot(depth[:, g3_id])
plt.show()
I found the example in (https://github.com/stoiver/anuga_warragamba_dam)
@chooron Your code looks fine. Though my guess is that you have not set the stage correctly to initially "fill" the reservoir. Also there could be a problem where the elevation DEM may create numerical dams which need to be filled before flowing into the actual reservoir.
@stoiver hello, I am back to working on this project, but I have found that when the flow rate is small, such as q=3e3, the water hardly flows forward. However, when I set the flow rate to be very large (q=3e5), the water flow can be observed, but this flow rate is too large and unrealistic. Additionally, I want the water to flow in a specific direction (since this is the process of opening the floodgates of a reservoir). Therefore, I would like to ask for your help again. How should I modify my code to address these issues?
@chooron happy to help. But am at present in the middle of a "hackathon" with the aim of implementing a GPU version of ANUGA. I'll should have some time to next week to have a quick look at your code. My guess is that you have numerical dams (due to mesh being too coarse in regions) which have to be filled before water flows downstream.