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)
14 Mar 23, 2022
36.9k Dec 28, 2022
3.1k Dec 31, 2022
1k Dec 29, 2022
349 Dec 25, 2022
8.1k Dec 30, 2022
0 Sep 05, 2021
1k Jan 08, 2023
1.6k Jan 04, 2023
6.4k Jan 02, 2023
16 Jun 30, 2022
424 Dec 15, 2022
137 Dec 22, 2022
180 Dec 28, 2022
9.9k Jan 08, 2023
915 Dec 29, 2022
1 Jan 03, 2022
133 Jan 04, 2023
62 Dec 08, 2022