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)
62 Dec 8, 2022
89 Jan 4, 2023
6 Nov 30, 2022
13 Sep 28, 2022
1 Jan 28, 2022
46 Dec 14, 2022
117 Dec 20, 2022
1 May 17, 2022
9 Aug 9, 2022
829 Jan 7, 2023
3 Mar 9, 2022
24 Nov 18, 2022
179 Dec 21, 2022
506 Dec 28, 2022
1.3k Jan 6, 2023
3 Feb 15, 2022
4 Nov 26, 2022
1 Dec 3, 2021
3 Sep 16, 2022
2 Sep 20, 2022