"""Example of shallow water wave equation analytical solution consists of a symmetrical converging channel with friction and bed slope. Specific methods pertaining to the 2D shallow water equation are imported from shallow_water for use with the generic finite volume framework Copyright 2004 Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray Geoscience Australia, ANU Specific methods pertaining to the 2D shallow water equation are imported from shallow_water for use with the generic finite volume framework Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the numerical vector named conserved_quantities. Existance of 'MacDonald_77541.tsh' assumed. """ #--------------- # Module imports import sys, string from os import sep sys.path.append('..'+sep+'pyvolution') from shallow_water import Transmissive_boundary, Reflective_boundary, \ Dirichlet_boundary from shallow_water import Domain, Add_value_to_region from pmesh2domain import pmesh_to_domain_instance from Numeric import zeros, Float #------- # Domain filename = 'MacDonald_91328_wall.tsh' print 'Creating domain from', filename domain = pmesh_to_domain_instance(filename, Domain) print 'Number of triangles = ', len(domain) #------------------------------------- # Provide file name for storing output domain.store = True #Store for visualisation purposes domain.format = 'sww' #Native netcdf visualisation format domain.set_name('MacDonald_first_order_wall') #-------------- # Bed Elevation fid = open('MacDonald_bed.xya') lines = fid.readlines() n_bed = len(lines) z_bed = zeros(n_bed, Float) x_bed = zeros(n_bed, Float) for line in range(n_bed): value = string.split(lines[line]) x_bed[line] = float(value[0]) z_bed[line] = float(value[1]) #------------------ # Set bed elevation def bed_elevation(x,y): n = x.shape[0] z = 0*x for i in range(n): # print i for j in range(n_bed-1): if x[i] >= x_bed[j] and x[i] < x_bed[j+1]: z[i] = z_bed[j] + (z_bed[j+1] - z_bed[j])/(x_bed[j+1] - x_bed[j])*(x[i] - x_bed[j]) return z # Bed elevation domain.set_quantity('elevation', bed_elevation) #----------------------------------------- # Set initial water surface stage which is # bed elevation plus an arbitrary 0.2 #--------------------------- # Add value to tagged region domain.set_region(Add_value_to_region('bed', 'stage', 0.2,location='unique vertices', initial_quantity='elevation')) #--------------------------- # Add value to tagged region domain.set_region(Add_value_to_region('floodplain', 'elevation', 2,location='unique vertices', initial_quantity='elevation')) #---------------------------------------------------------- # Decide which quantities are to be stored at each timestep domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] #------------------------------------------ # Reduction operation for get_vertex_values from anuga.pyvolution.util import mean domain.reduction = mean #-------------------- # Boundary Conditions upstream_depth = 5.035 downstream_depth = 1.5 discharge = 20 tags = {} tags['upstream'] = Dirichlet_boundary([upstream_depth, 2, 0.0]) tags['exterior'] = Reflective_boundary(domain) tags['downstream'] = Dirichlet_boundary([downstream_depth, 2, 0.0]) domain.set_boundary(tags) #--------- # Friction domain.set_quantity('friction', 0.02) # ------------------- # Set order of scheme domain.default_order = 1 domain.smooth = True #---------- # Evolution import time t0 = time.time() for t in domain.evolve(yieldstep = 100, finaltime = 1500.): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0)