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

Last change on this file since 2152 was 1354, checked in by chris, 20 years ago
File size: 3.0 KB
Line 
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   
7Specific methods pertaining to the 2D shallow water equation
8are imported from shallow_water
9for use with the generic finite volume framework
10
11Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
12numerical vector named conserved_quantities.
13
14Existence of file 'merimbula_interpolated.tsh' is assumed.
15"""
16
17#------------------------------
18# Setup Path and import modules
19import sys
20from os import sep, path
21sys.path.append('..'+sep+'pyvolution')
22
23from shallow_water import Domain, Reflective_boundary, File_boundary,\
24     Dirichlet_boundary, Wind_stress
25from pmesh2domain import pmesh_to_domain_instance
26from util import file_function, Polygon_function, read_polygon, inside_polygon
27from Numeric import zeros, Float, asarray
28from least_squares import Interpolation
29import time
30
31#-------
32# Domain
33filename = 'merimbula_interpolated_bathymetry.tsh'
34filename = 'merimbula_10834_bridge_refined_bathymetry.tsh'
35print 'Creating domain from', filename
36domain = pmesh_to_domain_instance(filename, Domain)
37print "Number of triangles = ", len(domain)
38
39#------------------------------------------
40# Reduction operation for get_vertex_values             
41from util import mean
42domain.reduction = mean
43
44domain.set_quantity('friction',0.03)
45#--------------------
46# Boundary conditions
47
48#---------------------------------------
49#   Tidal cycle recorded at Eden as open
50filename = 'Eden_tide_Sept03.dat'
51print 'Open sea boundary condition from ',filename
52Bf = File_boundary(filename, domain)
53
54#--------------------------------------
55#   All other boundaries are reflective
56Br = Reflective_boundary(domain)
57domain.set_boundary({'exterior': Br, 'open': Bf})
58
59#-----------
60# Wind field
61#   Format is time [DD/MM/YY hh:mm:ss], speed [m/s] direction (degrees)
62filename = 'Merimbula_Weather_data_Sept03_m_per_s.dat'
63print 'Wind field from ',filename
64F = file_function(filename, domain)
65domain.forcing_terms.append(Wind_stress(F))
66
67#--------------------------------
68# Initial water surface elevation
69domain.set_quantity('stage', -50.0)
70   
71#----------------------------------------------------------
72# Decide which quantities are to be stored at each timestep
73domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
74
75#-------------------------------------
76# Provide file name for storing output
77domain.store = True     #Store for visualisation purposes
78domain.format = 'sww'   #Native netcdf visualisation format
79domain.filename = 'Merimbula_2003_4days_dry'
80
81#----------------------
82# Set order of accuracy
83domain.default_order = 1
84domain.smooth = True
85print domain.set_inscribed_circle()
86
87#---------
88# Evolution
89t0 = time.time()
90yieldstep = 900
91finaltime = 588000*3
92for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
93    domain.write_time()
94   
95print 'That took %.2f seconds' %(time.time()-t0)   
Note: See TracBrowser for help on using the repository browser.