source: inundation/analytical solutions/Hobart.py @ 2152

Last change on this file since 2152 was 1367, checked in by chris, 20 years ago
File size: 3.1 KB
Line 
1"""Validation study of Merimbula lake using Pyvolution.
2Example of shallow water wave equation applied to
3Malpasset dam break simulation.
4
5   Copyright 2004
6   Christopher Zoppou, Stephen Roberts
7   Australian National University
8   
9Specific methods pertaining to the 2D shallow water equation
10are imported from shallow_water
11for use with the generic finite volume framework
12
13Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
14numerical vector named conserved_quantities.
15
16Existence of file 'Hobart_mesh.tsh' is assumed.
17"""
18
19###############################
20# Setup Path and import modules
21import sys
22from os import sep, path
23sys.path.append('..'+sep+'pyvolution')
24
25from shallow_water import Domain, Reflective_boundary, File_boundary,\
26     Dirichlet_boundary, Transmissive_boundary, Constant_height
27from pmesh2domain import pmesh_to_domain_instance
28from util import inside_polygon, read_polygon
29
30######################
31# Domain
32filename = 'Hobart_mesh.tsh'
33yieldstep = 1
34finaltime = 100
35   
36print 'Creating domain from', filename
37domain = pmesh_to_domain_instance(filename, Domain)
38print 'Number of triangles = ', len(domain)
39print 'Extent = ', domain.get_extent()
40
41
42domain.default_order = 1
43domain.smooth = True
44
45#domain.visualise = True
46
47
48#------------------------------
49# Boundary Conditions
50tags = {}
51tags['external'] = Reflective_boundary(domain) 
52domain.set_boundary(tags)
53
54#-----------------
55#Initial condition
56#   Throughout
57domain.set_quantity('stage', 0.)
58
59#    Within polygon region
60#        Set elevation to bed plus depth
61depth = 10.
62#        Read polygon boundary
63p0 = read_polygon('Hobart_recent_source_zone.xya')
64#        Get all centroid coordinates
65X = domain.get_centroid_coordinates()
66#        Find triangle indices which are within polygon boundary only
67indices = inside_polygon(X, p0)
68#        Get the bed elevation at the centroid of every triangle in polygon region
69
70zp = domain.get_quantity('elevation', location='centroids', indexes=indices)
71#        Set the stage to bed elevation plus depth for all triangle vertices in the polygon region
72
73domain.set_quantity('stage', zp+depth, location='centroids', indexes=indices) 
74
75
76#Alternative way which also works (I used it for testing - OMN)
77#
78
79#z = domain.get_quantity('elevation')
80#from copy import copy
81#w = copy(z)
82#for i in indices:
83#    w[i, :] = z[i, :] + depth
84#domain.set_quantity('stage', w)
85
86
87#-------------------------------------
88# Provide file name for storing output
89domain.store = True
90domain.format = 'sww'
91domain.filename = 'Hobart_first_order'
92
93#---------------------------------------------------------
94#Decide which quantities are to be stored at each timestep
95domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
96
97#----------------------------
98# Friction
99domain.set_quantity('friction', 0.05)
100         
101######################
102#Evolution
103import time
104t0 = time.time()
105for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
106    domain.write_time()
107   
108print 'That took %.2f seconds' %(time.time()-t0)
109
110   
Note: See TracBrowser for help on using the repository browser.