Hi, I am use ANUGA to simulate rainfall over DEM.
It works as I used uniform size mesh anuga.rectangular_cross_domain().
ncols = 648
nrows = 802
xllcorner = -52.17979305225
yllcorner = -1283978.9725515
cellsize = 290.34756669662
domain = anuga.rectangular_cross_domain(ncols, nrows, ncols*cellsize, nrows*cellsize, [xllcorner, yllcorner])
However, it gave warning as I used non-uniform mesh, and produce zero depth in the results, also delta t = 1000, the piece of code:
bounding_polygon = anuga.read_polygon('extent.csv') # Extent.
domain = anuga.create_domain_from_regions(bounding_polygon,
boundary_tags={'top': [0],
'right': [1],
'bottom': [2],
'left': [3]},
maximum_triangle_area=50000,
#interior_regions=[[bounding_polygon,500000]],
mesh_filename='tx_300m.msh',
use_cache=False,
verbose=True,
fail_if_polygons_outside=True)
The warning is:
/home/wlai/anaconda3/envs/anuga/lib/python2.7/site-packages/anuga/shallow_water/shallow_water_domain.py:2219: RuntimeWarning: invalid value encountered in less
negative_ids = num.where( num.logical_and((Stage.centroid_values - Elev.centroid_values) < 0.0 , tff > 0) )[0]
/home/wlai/anaconda3/envs/anuga/lib/python2.7/site-packages/anuga/file/sww.py:343: RuntimeWarning: invalid value encountered in greater_equal
storable_indices = num.array(w-z >= self.minimum_storable_height)
/home/wlai/anaconda3/envs/anuga/lib/python2.7/site-packages/anuga/file_conversion/sww2dem.py:197: RuntimeWarning: invalid value encountered in greater
q = fid.variables[name][:].flatten()
The whole code is also provided here. Thank you for your help.
"""Test with Tx300 DEM.
"""
#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
# Standard modules
import os
import time
import sys
import anuga
import numpy
from math import sin, pi, exp
time00 = time.time()
#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
anuga.asc2dem('tx_300m.asc')
anuga.dem2pts('tx_300m.dem')
# FIXME, tx300m Rainfall using non-uniform mesh doesnt work.
bounding_polygon = anuga.read_polygon('extent.csv') # Extent.
domain = anuga.create_domain_from_regions(bounding_polygon,
boundary_tags={'top': [0],
'right': [1],
'bottom': [2],
'left': [3]},
maximum_triangle_area=50000,
#interior_regions=[[bounding_polygon,500000]],
mesh_filename='tx_300m.msh',
use_cache=False,
verbose=True,
fail_if_polygons_outside=True)
'''
# Try uniform mesh. This works.
ncols = 648
nrows = 802
xllcorner = -52.17979305225
yllcorner = -1283978.9725515
cellsize = 290.34756669662
domain = anuga.rectangular_cross_domain(ncols, nrows, ncols*cellsize, nrows*cellsize, [xllcorner, yllcorner])
'''
# Print some stats about mesh and domain
print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
#print domain.statistics()
#------------------------------------------------------------------------------
# Setup parameters of computational domain
#------------------------------------------------------------------------------
domain.set_name('tx_300m') # Name of sww file
domain.set_datadir('.') # Store sww output here
domain.set_minimum_storable_height(0.01) # Store only depth > 1cm
domain.set_flow_algorithm('DE0')
#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
domain.set_quantity('friction', 0.03)
domain.set_quantity('elevation',
filename='tx_300m.pts',
use_cache=True,
verbose=True,
alpha=0.1)
h = 0.0
domain.set_quantity('stage', expression='elevation + %f' % h) ## FIXME, set depth?
#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
bc_wall = anuga.Reflective_boundary(domain) # Solid reflective wall
bc_out = anuga.Transmissive_boundary(domain)
# Associate boundary tags with boundary objects
domain.set_boundary({'left': bc_out, 'right': bc_out, 'top': bc_out, 'bottom': bc_out})
#---------------------------------------------------------------------------------------
# Rainfall WL. Note rainfall file is daily accumulated rainfall in inch/day.
anuga.asc2dem('tx_rainfall.asc')
anuga.dem2pts('tx_rainfall.dem')
ncols = 41
nrows = 50
xllcorner = -4757.5382496486
yllcorner = -1288482.8472377
cellsize = 4763
rainfall_domain = anuga.rectangular_cross_domain(ncols, nrows, ncols*cellsize, nrows*cellsize, [xllcorner, yllcorner])
# Print some stats about mesh and domain
#print 'Number of triangles = ', len(rainfall_domain)
#print 'The extent is ', rainfall_domain.get_extent()
#print rainfall_domain.statistics()
rainfall_domain.set_quantity('stage',
filename='tx_rainfall.pts',
use_cache=True,
verbose=True,
alpha=0.1)
#print 'rainfal domain ', dir(rainfall_domain)
conv = 2.9398e-7 # convert inch/day to meter/sec.
inter_pts = domain.get_vertex_coordinates()
rf_accu = rainfall_domain.get_quantity('stage').get_values(interpolation_points=inter_pts)
rf_accu = numpy.reshape(rf_accu, (domain.get_number_of_triangles(), 3))
rf_rate = conv * rf_accu
print 'number of inter_pts', len(inter_pts), ',size of inter_pts', sys.getsizeof(inter_pts), ',shape:', inter_pts.shape
print 'number of rf_rate', len(rf_rate), ',size of rf_rate', sys.getsizeof(rf_rate), ',shape:', rf_rate.shape
#print 'domain.eleveation shape : ', len(domain.quantities['elevation']), ', col', len((domain.quantities['elevation'])[0])
'''
print domain.get_number_of_nodes()
print domain.get_number_of_triangles()
print domain.get_unique_vertices()
print rainfall_domain.get_number_of_nodes()
print rainfall_domain.get_number_of_triangles()
'''
#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
dt = 60 # Timestep in seconds.
rf_rate = abs(rf_rate * dt)
for t in domain.evolve(yieldstep=dt, finaltime=1200.0):
print domain.timestepping_statistics()
domain.add_quantity('stage', rf_rate)
# Print time.
print 'That took %.2f seconds' %(time.time()-time00)
# Create inundation map.
anuga.sww2dem('tx_300m.sww', 'inundation.asc', quantity='stage-elevation', cellsize=300,reduction=max,verbose=True)
15 May 07, 2022
2 Nov 14, 2021
6.3k Jan 03, 2023
276 Dec 19, 2022
15.7k Jan 07, 2023
186 Dec 23, 2022
62 Dec 08, 2022
1k Dec 29, 2022
3.6k Dec 27, 2022
424 Dec 15, 2022
180 Dec 28, 2022
535 Jan 04, 2023
0 Sep 05, 2021
613 Dec 16, 2022
944 Dec 28, 2022
3.2k Jan 08, 2023
1.6k Jan 04, 2023
32 Dec 14, 2022
1.8k Dec 31, 2022