source: production/gippsland_2005/run_gippsland2.py @ 2149

Last change on this file since 2149 was 2149, checked in by duncan, 18 years ago

'scenario'

File size: 4.6 KB
Line 
1
2#Convention for strings representing files
3# #_file has the extention
4#           #name does not have the extension
5
6import time
7
8from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
9     dem2pts, asc_csiro2sww
10from pyvolution.pmesh2domain import pmesh_to_domain_instance
11from caching import cache
12from create_mesh2 import create_mesh
13from pyvolution.shallow_water import Domain, Reflective_boundary,\
14     File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary
15from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA
16
17# Import this last!  It stuffs up the loading of c extensions otherwise.
18import project
19
20#Data preparation
21#Convert ASC 2 DEM 2 PTS using source data and store result in source data
22# just interested in the boundary/ mesh intereaction first off.
23#cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
24#      dependencies = [demname + '.asc'],
25#      verbose = True)
26#      #evaluate = True)
27
28cache(dem2pts, project.demname[:-4], {'verbose': True},
29      dependencies = [project.demname]     
30      ,verbose = False
31      )
32
33
34# Convert CSIRO boundary
35# NOTE: This function is dependent on a lot of files, and I haven't listed
36# them, since there is so many and I don't think they will be changing.
37cache(asc_csiro2sww,
38      (project.bath_dir,
39       project.elevation_dir,
40       project.ucur_dir,
41       project.vcur_dir,
42       project.boundaryname),
43       {'verbose': True,
44       'minlat': -38.25,
45       'maxlat': -37.7,
46       'minlon': 147.0,
47       'maxlon': 148.25,
48       'mean_stage': project.tide,
49       'zscale': 1,                 #Enhance storm surge
50       'fail_on_NaN': False}
51      #,evaluate = True                 
52      ,verbose = False
53      )
54
55meshname = project.meshname
56pointname = project.pointname
57mesh_elevname = project.mesh_elevname
58outputname = project.outputname
59
60t0 = time.time()
61#Create the mesh without elevation data
62# 100000 very course - 16,827 triangle mesh
63
64meshname, triagle_count = cache(create_mesh,(1000000),
65                                {'factor':1000,
66                                 'mesh_file':meshname,
67                                 'triangles_in_name':True}
68                                ,dependencies = ['create_mesh2.py']
69                                #,evaluate = True     
70                                )
71
72mesh_elevname = mesh_elevname[:-9] + '_' + str(triagle_count) +  \
73                mesh_elevname[-9:]
74outputname = outputname[:-4] + '_' + str(triagle_count) + outputname[-4:]
75
76
77#Add elevation data to the mesh
78#fit_to_mesh_file(mesh_file, point_file, mesh_output_file, DEFAULT_ALPHA,
79#                 verbose=True, expand_search=True, precrop=True)     
80cache(fit_to_mesh_file,(meshname,
81                 pointname,
82                 mesh_elevname,
83                 DEFAULT_ALPHA),
84      {'verbose': True,
85       'expand_search': True,
86       'precrop': True}
87      ,dependencies = [meshname, pointname]
88      #,evaluate = True     
89      ,verbose = False
90      )
91 
92print 'Initialising the mesh took %.2f seconds' %(time.time()-t0) 
93
94#Setup domain
95domain = cache(pmesh_to_domain_instance, (mesh_elevname, Domain),
96               dependencies = [mesh_elevname]                   
97               ,verbose = False
98               )               
99
100
101domain.set_name(project.basename + '_%d' %triagle_count)
102domain.set_datadir(project.outputdir)
103domain.store = True
104domain.quantities_to_be_stored = ['stage']
105
106print 'Number of triangles = ', len(domain)
107print 'The extent is ', domain.get_extent()
108
109#Setup Initial Conditions
110domain.set_quantity('friction', 0)
111domain.set_quantity('stage', project.tide)
112
113
114#Setup Boundary Conditions
115print domain.get_boundary_tags()
116
117print "****** run  ****"
118print "project.boundaryname",project.boundaryname
119print "**********"
120
121Bf = File_boundary(project.boundaryname, domain, verbose = True)
122#domain.starttime = 3000  #Obtained from MOST
123domain.starttime = 0  #Obtained from MOST
124
125Br = Reflective_boundary(domain)
126Bt = Transmissive_boundary(domain)
127Bd = Dirichlet_boundary([project.tide,0,0])
128#Bw = Time_boundary(domain=domain,
129#                   f=lambda t: [(60<t<660)*4, 0, 0])
130
131domain.set_boundary( {'back': Br,'side': Bf, 'ocean': Bf} ) #CSIRO storm surge
132
133#Evolve
134t0 = time.time()
135
136#for t in domain.evolve(yieldstep = 60, finaltime = 30000):
137for t in domain.evolve(yieldstep = 10, finaltime = 30):
138#                       skip_initial_step = True):
139    domain.write_time()
140    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')       
141
142print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.