Drainage or Velocity Issue

[tyvaughn]

Hi all--

I am running the below code to simulate flooding in a 49km^2 watershed with a very shallow slope. Although I am simulating a pretty large event (30cm in 8 hours) the water levels that are returned to me are half as large as they should be. Does anyone have any ideas? Thanks in advance.

############################################################################
# Import ANUGA
############################################################################
import sys
default_stdout = sys.stdout
default_stderr = sys.stderr

import anuga
reload(sys)
sys.stdout = default_stdout
sys.stderr = default_stderr

from netCDF4 import Dataset

############################################################################
# Convert elevation raster from ASCII to PTS file
############################################################################

filename_root = 'data/SCDEM'

anuga.asc2dem(filename_root + '.asc', use_cache = False, verbose = True)
anuga.dem2pts(filename_root + '.dem', use_cache = False, verbose = True)

############################################################################
# Import boundary polygon (CSV of node UTM coordinates) and set boundary tags
############################################################################

bounding_polygon = anuga.read_polygon('data/extent.csv')
boundary_tags = {'all':[0,1,2,3,4,5,6,7,8]} 

############################################################################
# Create the computational domain
############################################################################

domain = anuga.create_domain_from_regions(bounding_polygon = bounding_polygon,
                                          boundary_tags = boundary_tags,
                                          maximum_triangle_area = 1600,
                                          mesh_filename = filename_root + '.msh')

############################################################################
# Set initial conditions
############################################################################

domain.set_quantity('elevation', 
                    filename = filename_root + '.pts')

domain.set_quantity('stage', expression='elevation')   # Dry initial condition

domain.set_quantity('friction', 0.035)

############################################################################
# Set other domain options
############################################################################

domain.set_name('Practice')         # Name of sww file
swwFile = 'Practice.sww'

domain.set_quantities_to_be_stored({'elevation': 1,
                                    'stage': 2,
                                    'xmomentum': 2,
                                    'ymomentum': 2})

############################################################################
# Set boundary conditions
############################################################################

min_elev = domain.quantities['elevation'].vertex_values.min()

Bo = anuga.Dirichlet_boundary([min_elev, 0., 0.]) # outlet

domain.set_boundary({'all':Bo})

############################################################################
# Rainfall 
############################################################################

from anuga.operators.rate_operators import Rate_operator

initial_rate = 0.0000103 #3.5cm/hr

rain_operator = Rate_operator(domain, rate = initial_rate)

############################################################################
# Evolve the domain
############################################################################

timestep = 150

for t in domain.evolve(yieldstep = timestep, finaltime = 3600*6):
    print domain.timestepping_statistics()

    if t > 3600:
        rain_operator.set_rate(rate = 0.0)

############################################################################
# FID
############################################################################
fid = Dataset(swwFile, mode='r')

x = fid.variables['x'][:]
y = fid.variables['y'][:]

elev = fid.variables['elevation'][:]
stage = fid.variables['stage'][:]
depth = stage - elev

xmom = fid.variables['xmomentum'][:]
xvel = (xmom * depth) / (depth**2 + 0.0001)

ymom = fid.variables['ymomentum'][:]
yvel = (ymom * depth) / (depth**2 + 0.0001)

fid.close()

############################################################################
# Create Time Series
############################################################################

gauge_filename = 'gauges/gauges.csv'

try:
    anuga.sww2csv_gauges((swwFile),
    gauge_filename,
    quantities=['depth'],
    verbose=True)

except:
    print 'Failed to process gauges'

[rudyvandrie]

The rate for rainfall in ANUGA is mm/sec and with a factor of 10^-3 .
So 0.00000103
You want 35mm/hr = 35/3600 = 0.0097mm/sec