1 | """Validation study of Merimbula lake using Pyvolution. |
---|
2 | |
---|
3 | Copyright 2004 |
---|
4 | Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray |
---|
5 | Geoscience Australia, ANU |
---|
6 | |
---|
7 | Specific methods pertaining to the 2D shallow water equation |
---|
8 | are imported from shallow_water |
---|
9 | for use with the generic finite volume framework |
---|
10 | |
---|
11 | Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the |
---|
12 | numerical vector named conserved_quantities. |
---|
13 | |
---|
14 | Existence of file 'merimbula_10785_1.tsh' is assumed. |
---|
15 | """ |
---|
16 | |
---|
17 | #------------------------------ |
---|
18 | # Setup Path and import modules |
---|
19 | import sys |
---|
20 | from os import sep, path |
---|
21 | sys.path.append('..'+sep+'pyvolution') |
---|
22 | |
---|
23 | from shallow_water import Domain, Reflective_boundary, File_boundary,\ |
---|
24 | Dirichlet_boundary, Wind_stress |
---|
25 | from pmesh2domain import pmesh_to_domain_instance |
---|
26 | from util import file_function, Polygon_function, read_polygon, inside_polygon |
---|
27 | from Numeric import zeros, Float, asarray |
---|
28 | from least_squares import Interpolation |
---|
29 | import time |
---|
30 | |
---|
31 | |
---|
32 | #------- |
---|
33 | # Domain |
---|
34 | filename = 'merimbula_10785_1.tsh' |
---|
35 | print 'Creating domain from', filename |
---|
36 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
37 | print "Number of triangles = ", len(domain) |
---|
38 | |
---|
39 | #------------------------------------------ |
---|
40 | # Reduction operation for get_vertex_values |
---|
41 | from util import mean |
---|
42 | domain.reduction = mean |
---|
43 | |
---|
44 | |
---|
45 | domain.set_quantity('friction',0.03) |
---|
46 | #-------------------- |
---|
47 | # Boundary conditions |
---|
48 | |
---|
49 | #--------------------------------------- |
---|
50 | # Tidal cycle recorded at Eden as open |
---|
51 | filename = 'Eden_tide_Sept03.tms' |
---|
52 | print 'Open sea boundary condition from ',filename |
---|
53 | Bf = File_boundary(filename, domain) |
---|
54 | |
---|
55 | #-------------------------------------- |
---|
56 | # All other boundaries are reflective |
---|
57 | Br = Reflective_boundary(domain) |
---|
58 | domain.set_boundary({'exterior': Br, 'open': Bf}) |
---|
59 | |
---|
60 | #----------- |
---|
61 | # Wind field |
---|
62 | # Format is time [DD/MM/YY hh:mm:ss], speed [m/s] direction (degrees) |
---|
63 | #filename = 'Merimbula_Weather_data_Sept03_m_per_s.dat' |
---|
64 | #print 'Wind field from ',filename |
---|
65 | #F = file_function(filename, domain) |
---|
66 | #domain.forcing_terms.append(Wind_stress(F)) |
---|
67 | |
---|
68 | #-------------------------------- |
---|
69 | # Initial water surface elevation |
---|
70 | domain.set_quantity('stage', -50.0) |
---|
71 | |
---|
72 | #---------------------------------------------------------- |
---|
73 | # Decide which quantities are to be stored at each timestep |
---|
74 | domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] |
---|
75 | |
---|
76 | #------------------------------------- |
---|
77 | # Provide file name for storing output |
---|
78 | domain.store = True #Store for visualisation purposes |
---|
79 | domain.format = 'sww' #Native netcdf visualisation format |
---|
80 | from normalDate import ND |
---|
81 | domain.filename = 'Merimbula_2003_4days_dry_%s'%ND() |
---|
82 | |
---|
83 | print domain.filename |
---|
84 | |
---|
85 | #---------------------- |
---|
86 | # Set order of accuracy |
---|
87 | domain.default_order = 1 |
---|
88 | domain.smooth = True |
---|
89 | print domain.use_inscribed_circle |
---|
90 | |
---|
91 | #--------- |
---|
92 | # Evolution |
---|
93 | t0 = time.time() |
---|
94 | domain.visualise = False |
---|
95 | yieldstep = 60 |
---|
96 | finaltime = 3600*24*3 |
---|
97 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
98 | domain.write_time() |
---|
99 | |
---|
100 | print 'That took %.2f seconds' %(time.time()-t0) |
---|