"""Validation study of Merimbula lake using Pyvolution. Example of shallow water wave equation applied to Malpasset dam break simulation. Copyright 2004 Christopher Zoppou, Stephen Roberts Australian National University 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. Existence of file 'Hobart_mesh.tsh' is assumed. """ ############################### # Setup Path and import modules import sys from os import sep, path sys.path.append('..'+sep+'pyvolution') import Numeric from shallow_water import Domain, Reflective_boundary, File_boundary,\ Dirichlet_boundary, Transmissive_boundary from pmesh2domain import pmesh_to_domain_instance from util import Polygon_function, read_polygon ###################### # Domain filename = 'hobart_mesh.tsh' filename = 'building_boundary_Hobart_58868_triangles12May2005_bathymetry_0.1.tsh' yieldstep = 0.1 finaltime = 2000 print 'Creating domain from', filename domain = pmesh_to_domain_instance(filename, Domain) print "Number of triangles = ", len(domain) domain.default_order = 1 domain.smooth = True # Turn on visualiser domain.initialise_visualiser() # Colours for Hobart domain.visualiser.stage_color = (0.0,0.38,0.1) domain.visualiser.friction_color = (0.1,0.99,0.1) domain.visualiser.other_color = (0.99,0.4,0.1) #------------------------------ # Boundary Conditions tags = {} tags['exterior'] = Reflective_boundary(domain) tags['Building'] = None domain.set_boundary(tags) #--------------------------------------------------------- #Decide which quantities are to be stored at each timestep domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] #----------------- #Initial condition p0 = read_polygon('Hobart_recent_source_zone.xya') domain.set_quantity('stage',Polygon_function([(p0,10.0)])) stage_vv = domain.quantities['stage'].vertex_values elevation_vv = domain.quantities['elevation'].vertex_values new_stage_vv = stage_vv + elevation_vv domain.set_quantity('stage', new_stage_vv) #------------------------------------- # Provide file name for storing output domain.store = True domain.format = 'sww' domain.filename = 'Hobart_first_order' #---------------------------- # Friction domain.set_quantity('friction', 0.033) Friction = domain.quantities['friction'] for triangle in domain.tagged_elements['Building']: #print 'key %s triangle %g \n'%(key, triangle) Friction.centroid_values[triangle] = 100.0 Friction.extrapolate_first_order() ###################### #Evolution import time t0 = time.time() for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): domain.write_time() domain.visualiser.update_quantity_color('friction',(1.0,0.0,0.0)) print 'That took %.2f seconds' %(time.time()-t0)