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)
1 Dec 16, 2021
2 Nov 14, 2021
915 Dec 29, 2022
48.9k Jan 03, 2023
133 Jan 04, 2023
35 Dec 13, 2022
535 Jan 04, 2023
1.6k Jan 04, 2023
443 Jan 01, 2023
3.9k Jan 05, 2023
3.6k Dec 27, 2022
267 Jan 01, 2023
6.4k Jan 02, 2023
1k Jan 08, 2023
9.9k Jan 08, 2023
349 Dec 25, 2022
705 Jan 05, 2023
15 May 07, 2022
8.1k Dec 30, 2022
37 Dec 21, 2022