source: anuga_work/development/cairns_2006/run_cairns.py @ 4113

Last change on this file since 4113 was 3944, checked in by sexton, 18 years ago

(i) update sww2timeseries so can handle gauges which don't fall within sww domain (ii) script for generating report for dampier based on sww files from parallel setup

File size: 4.1 KB
RevLine 
[2229]1"""
2Author: John Jakeman
3Created: 12/12/2005
4Main Cairns script using new interface
5"""
6
7#-------------------------------
8# Module imports
9#-------------------------------
10import sys, os
[3944]11from anuga.shallow_water import Domain, Reflective_boundary,\
[2312]12     File_boundary, Transmissive_Momentum_Set_Stage_boundary,\
13     Transmissive_boundary, Dirichlet_boundary
[3944]14#from anuga.pyvolution.mesh_factory import rectangular_cross
15#from anuga.pmesh2domain import pmesh_to_domain_instance
16from anuga.pmesh.mesh_interface import create_mesh_from_regions
[2229]17from Numeric import array, zeros, Float, allclose
[3434]18import cairns_project as project
19#import project
[2229]20from caching import cache
[3514]21from anuga.utilities.polygon import read_polygon, Polygon_function
[2229]22
[2312]23from conversion import convert_latlon_to_xycoords
24
[2229]25#-------------------------------
26# Domain
27#-------------------------------
28print 'Creating domain'
[2312]29
[2270]30"""
31M = 36 #float(len1)/m
32N = 10 #float(len2)/n
[2229]33
34points, elements, boundary = rectangular_cross(M, N, len1=6393693.9504678631,
35                                               len2=1781111.2941870166,
36                                               origin = (0.0, 0.0))
37
38#points, elements, boundary = rectangular_cross(N, M, len1=58.0, len2=1.0, origin=(145.0,-24.0))
39                                               
[2270]40domain = cache(Domain, (points, elements, boundary))
[2229]41
[2270]42"""
[2312]43
[3944]44#cache(pmesh_to_domain_instance,
45#      (project.mesh_filename, Domain),
46#      dependencies = [project.mesh_filename])
47domain = Domain(project.mesh_filename, use_cache = True, verbose = True)
[2270]48
[2229]49domain.check_integrity()
50print 'Number of triangles = ', len(domain)
51print 'The extent is ', domain.get_extent()
52
53#---------------------------------------------------------
54#Decide which quantities are to be stored at each timestep
55domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
56
57#-------------------------------
58# Initial Conditions
59#-------------------------------
60print 'Initial values'
61
[2312]62latlong_origin = 145, -24
63r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin)
64r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin)
65r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin)
66r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin)
67
68#p0 = read_polygon('cairns_fault.xya')
[2229]69#p0 = read_polygon('cairns_fault_degrees.xya')
[2312]70p0 = [r0, r1, r2, r3]
71domain.set_quantity('stage',Polygon_function([(p0,1.0)]))
[2229]72
73#print domain.quantities['stage'].vertex_values
74
75domain.set_quantity('elevation',
[2270]76                    filename = project.bathymetry_filename[:-4] + '.xya',
[2229]77                    alpha = 10.0,
78                    verbose = True,
79                    use_cache = True)
80
81# Need to increase elevation by the same amount as the stage
82# At present unsure how to achieve this
83#domain.set_quantity('elevation',Polygon_function([domain.quantities['elevation'].vertex_values+1000.0]))
84
85domain.set_quantity('friction', 0.0)
86
87#-------------------------------
88# Boundary conditions
89#-------------------------------
90print 'Boundaries'
91
92#Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function)
[2312]93#Bd = Dirichlet_boundary([0,0,0])
[2270]94Br = Reflective_boundary(domain)
95domain.set_boundary({'border': Br, 'border': Br, 'border': Br, 'border': Br})
[2229]96
[2270]97#Bt = Transmissive_boundary(domain)
[2312]98#domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
[2229]99
100#-------------------------------
101# Setup domain runtime parameters
102#-------------------------------
[2270]103
[2312]104#domain.visualise = True
105#domain.visualise_color_stage = True
[2229]106
[2312]107#Set domain.visualise = False for large problem sizes
108domain.visualise = False
109domain.visualise_color_stage = False
[2270]110
[2229]111base = os.path.basename(sys.argv[0])
[3846]112basename, _ = os.path.splitext(base)
113domain.set_name(basename)
[2229]114domain.default_order = 2
115domain.store = True    #Store for visualisation purposes
116
117#-------------------------------
118# Evolve
119#-------------------------------
120import time
121t0 = time.time()
[2270]122yieldstep = 60.0
[2312]123finaltime = 21600.0
[2229]124
125for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
126    domain.write_time()
127
128print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.