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)
349 Dec 25, 2022
37 Dec 21, 2022
443 Jan 01, 2023
705 Jan 05, 2023
3.6k Dec 27, 2022
15.7k Jan 07, 2023
1.8k Dec 31, 2022
15 May 07, 2022
7.9k Jan 01, 2023
6.4k Jan 02, 2023
424 Dec 15, 2022
1k Dec 29, 2022
944 Dec 28, 2022
1.1k Jan 06, 2023
36.9k Dec 28, 2022
3 Aug 04, 2021
276 Dec 19, 2022
0 Sep 05, 2021
186 Dec 23, 2022