source: anuga_work/development/ted_rigby/model_run.py @ 7475

Last change on this file since 7475 was 5997, checked in by ole, 15 years ago

Cosmetics

File size: 21.3 KB
Line 
1"""
2----------------------------------------------------------------- model_run.py
3ANUGA RURAL FLOOD MODELLING TEMPLATE (Domain Construction & Execution)
4Script for building the model domain and setting initial conditions and
5(temporally varying)  boundary conditions appropriate to the particular
6scenario-event before running (evolving) the simulation.
7
8Most model data for the run is read in or otherwise referenced in
9model_data.py and imported into model_run.py to create the model mesh and
10domain before running the simulation.
11
12Code Version : 1.00 June 2008 Initial release
13               1.01 December 2008 Vble_n recoded to speed execution
14Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au
15
16-------------------------------------------------------------------------------
17"""
18
19#------------------------------------------------------------------------------
20# Import necessary modules
21#------------------------------------------------------------------------------
22
23# Standard python modules
24import os
25import time
26import sys
27from Numeric import zeros, Float
28
29# Related ANUGA modules
30from anuga.shallow_water import Domain
31from anuga.shallow_water import Reflective_boundary
32from anuga.shallow_water import Dirichlet_boundary
33from anuga.shallow_water import Time_boundary
34from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
35from anuga.shallow_water.shallow_water_domain import Inflow
36from anuga.abstract_2d_finite_volumes.util import file_function
37from anuga.abstract_2d_finite_volumes.quantity import Quantity
38from anuga.pmesh.mesh_interface import create_mesh_from_regions
39from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function, inside_polygon
40from anuga.alpha_shape.alpha_shape import alpha_shape_via_files
41
42
43# Model specific imports
44import model_data                # Defines and provides the scenario-event specific data
45import get_tuflow_data           # Provides the routines to read in tuflow specific data
46
47#####################################################################################
48#                                                                                   #
49#                     S I M U L A T I O N    M A I N L I N E                        #
50#                                                                                   #
51#####################################################################################
52
53t0 = time.time()                 # record start time of this run (secs)
54print model_data.basename+' >>>> Comencing model_run.py code execution at t= %.2f hours'  %(float(time.time()-t0)/3600.0)
55
56#####################################################################################
57#                                                                                   #
58#                          CREATE THE TRIANGULAR MESH                               #
59#                                                                                   #
60#####################################################################################
61
62# Note: The mesh is created based on overall clipping polygon with a tagged 
63# boundary and interior regions as defined in model_data.py along with       
64# resolutions (maximal area per triangle) for each interior polygon       
65
66print model_data.basename+' >>>> Creating variable mesh within bounding_poly reflecting regional resolutions'
67# Create the mesh reflecting the bounding (default) and internal region mesh resolutions set in model_data.py
68# with boundary tags as defined in model_data.py
69
70create_mesh_from_regions(model_data.bounding_polygon,
71                         boundary_tags = model_data.boundary_strings,
72                         maximum_triangle_area = model_data.default_res,
73                         minimum_triangle_angle=30.0,
74                         interior_regions = model_data.interior_resregions,
75                         filename = model_data.mesh_filename,
76                         use_cache=True,
77                         verbose=model_data.anuga_verbose)
78
79
80print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
81
82######################################################################################
83#                                                                                    #
84#                       CREATE THE COMPUTATIONAL DOMAIN                              #
85#                                                                                    #
86######################################################################################
87
88print model_data.basename+' >>>> Creating domain from mesh and setting domain parameters'
89# create domain from mesh meshname
90domain = Domain(model_data.mesh_filename, use_cache=True, verbose=model_data.anuga_verbose) 
91print domain.statistics()                                # confirm what has been done
92
93# set base domain parameters
94domain.set_name(model_data.basename)                     # already created in model_data.py
95domain.set_datadir(model_data.results_dir)                 # already created and/or checked in model_data.py
96domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
97domain.set_minimum_storable_height(0.01)                 # Remove very shallow water depths from sww
98domain.set_maximum_allowed_speed(8.0)                    # Allow Maximum velocity of 8m/s....   
99
100######################################################################################
101#                                                                                    #
102#                       ASSIGN THE DOMAIN INITIAL CONDITIONS                         #
103#                                                                                    #
104######################################################################################
105
106#-------------------------------------------------------------------------------------
107# This section assigns the initial (t=0) conditions for the domain
108# Domain initial conditions can be fixed, read from a file, or from a function
109# Note also that anuga permits almost all initial conditions to subsequently vary
110# with time, so while the following conditions apply to the whole domain at t=0
111# all can be modified within the evolve time loop. In this way scour, deposition,
112# blockage, changes in roughness etc with time can be explicitly included in a simulation
113#-------------------------------------------------------------------------------------
114
115print model_data.basename+' >>>> Assigning initial conditions to %s domain' %(model_data.basename)
116
117# Elevation data may be assigned to the domain from a csv or (NetCDF) pts file.
118
119# For a smaller csv that will be read in repeatedly, best to use geospatial_data to readin and
120# export_point_file to save the points permanently in a NetCDF *.pts file. This runs faster
121# than the block read of a csv - but - is limited to about 6 million points.
122# Alternatively - by directly assigning elevation from a csv file can take advantage of the block
123# read in domain.set_quantity so the elevation dataset can be virtually unlimited in size.
124
125# In this simulation the size of the elevation dataset is to large for conversion to a (NetCDF)
126# pts format so data is read in from the raw csv file.
127
128domain.set_quantity('elevation', 
129                    filename = model_data.elev_filename,
130                    use_cache=True,
131                    verbose=model_data.anuga_verbose,
132                    alpha=0.01)
133
134
135# Set the intial water level (stage) across the domain             
136InitialWaterLevel = model_data.InitialWaterLevel             # At t=0 water surface is in this model at initial lake level (0.3m)
137if domain.get_quantity('elevation') < InitialWaterLevel:     # This can create ponds if land is below initial tailwater level
138    domain.set_quantity('stage', InitialWaterLevel)          # set initial stage at initial tailwater level (wet)
139else:       
140    domain.set_quantity('stage', expression='elevation')     # else set initial stage at land level (dry)
141
142
143# Assign uniform initial friction value to domain mesh centroids as initial condition
144# Note: Friction will be re-assessed  based on roughness data read in from tuflow
145# 2d_mat files and computed depth during evolve time loop!!!!
146domain.set_quantity('friction', model_data.default_n, location = 'centroids')
147
148
149
150
151#######################################################################################
152#                                                                                     #
153#             ASSIGN THE (TEMPORAL) SUBAREAL INFLOWS  TO THE DOMAIN                   #
154#                                                                                     #
155#######################################################################################
156
157#--------------------------------------------------------------------------------------
158# This section reads the local or total temporal inflows read in in model_data.py from
159# the Tuflow ts1 files and stored as NETCDF tms file in the tms_files subdirectory,
160# into the domain.
161#
162# Note: The current ANUGA inflow function 'pours' flow onto the domain within a circle
163# or polygonal area specified by the user at the specified temporal rate.
164# EHR -- the above inflows should eventually be augmented by  inflows at the boundary
165# (total as distinct from local hydrographs) to better correspond with conventional
166# flood modelling.
167#--------------------------------------------------------------------------------------
168
169print model_data.basename+" >>>> Assigning the inflow hydrographs from the (%s) tms files directory" %(model_data.tms_dir)
170
171
172inflow_fields=[]
173# read in each line of inflow_hydrographs[tms_filename,xcoord,ycoord] and assign as inflow hydrograph to domain
174for i in range(len(model_data.inflow_hydrographs)) :
175   
176    tms_filename = model_data.inflow_hydrographs[i][0]                   # extract the inflow hydro tms_filename
177    xcoord = model_data.inflow_hydrographs[i][1]                         # extract the inflow hydro xcoord
178    ycoord = model_data.inflow_hydrographs[i][2]                         # extract the inflow hydro ycoord
179    print model_data.basename+' >>>> Reading inflow hydrograph from %s at %7.2f %7.2f and appending to forcing terms ' %(tms_filename,xcoord,ycoord)
180       
181    # convert time series to temporal function
182    flowrate = file_function(tms_filename, 
183                             quantities='value',
184                             boundary_polygon=None)
185    hydrograph = Inflow(domain,
186                        center=(xcoord,ycoord),
187                        radius=10, 
188                        rate=flowrate)
189
190    # append hydrograph to domain
191    domain.forcing_terms.append(hydrograph)
192
193print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.inflow_hydrographs),float(time.time()-t0)/3600.0)
194
195
196#######################################################################################
197#                                                                                     #
198#        ASSIGN THE (TEMPORAL) DOWNSTREAM BOUNDARY CONDITIONS TO THE DOMAIN           #
199#                                                                                     #
200#######################################################################################         
201
202
203print model_data.basename+' >>>> Assigning the DSBC from the (%s) tms files directory ' %(model_data.tms_dir)
204print model_data.basename+' >>>> Available boundary tags are ', domain.get_boundary_tags()
205
206# convert dsbc time series to function of time
207tms_filename=model_data.dsbc_hydrographs[0]
208lake_stage  = file_function(tms_filename,quantities='value',boundary_polygon=None)
209# Note: will need to be modified if more than one dsbc is to be applied
210print model_data.basename + ' >>>> Applying dsbc %s to model at the Lake_bdry ' %(tms_filename)
211
212MaMt_bdry   = Reflective_boundary(domain)    # all reflective except lake as using inflow()
213Mriv_bdry   = Reflective_boundary(domain)    # Really redundant untill get boundary inflows
214Cbck_bdry   = Reflective_boundary(domain) 
215WFra_bdry   = Reflective_boundary(domain) 
216EFra_bdry   = Reflective_boundary(domain) 
217
218lake_bdry   = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=lake_stage)
219closed_bdry = Reflective_boundary(domain)                         
220
221# boundary conditions for current scenario
222domain.set_boundary({'Lake': lake_bdry,
223                     'MaMt': MaMt_bdry,
224                     'Mriv': Mriv_bdry,
225                     'Cbck': Cbck_bdry,
226                     'WFra': WFra_bdry,
227                     'EFra': EFra_bdry,
228                     'exterior': closed_bdry })  # this relates to the residue of the bounding poly not explicitly asigned above
229   
230print model_data.basename+' >>>> Completed assignment of %i dsbc boundary conditions to domain at t= %.2f hours ' %(len(model_data.dsbc_hydrographs),float(time.time()-t0)/3600.0)
231
232
233#######################################################################################
234#                                                                                     #
235#  ASSIGN THE (TEMPORALLY DEPTH DEPENDENT) SURFACE ROUGHNESS (FRICTION) PARAMETERS    #
236#                                                                                     #
237#######################################################################################
238
239#---------------------------------------------------------------------------------------
240# This section creates an array of surface roughness data for each element in the domain
241# ready for later use in the evolve loop where the roughness data is  used with depth to
242# compute a depth dependent friction at each yieldstep (output timestep to the swww file).
243#
244# Note that because the polys read in should but may not cover the entire bounding poly area
245# It is prudent to initially declare a default material (as in Tuflow) applying the default
246# to the whole domain before applying (overwriting) with the spatial roughness data from the
247# polys obtained from the 2d_mat mid/mif files!!!!
248#--------------------------------------------------------------------------------------
249
250print model_data.basename+' >>>> Creating and assigning surface roughness data for use in evolve'
251
252N = len(domain)
253# Create the new roughness array to hold the elemental roughness data
254# In this array of Nx5 values
255#       n0 -- is the equivalent fixed roughness value 
256#       d1 -- is the depth below which n1 applies
257#       n1 -- is the roughness applying below d1
258#       d2 -- is the depth above which n2 applies
259#       n2 -- is the roughness applying above d2 
260
261# Create and assign zeros to the roughness array
262material_variables = zeros((N, 5), Float)               
263
264# Set default material variables n0, d1, n1, d2, n2 in the roughness array
265default_n0 = 0.040
266default_d1 = 0.3
267default_n1 = 0.050
268default_d2 = 1.5
269default_n2 = 0.030
270
271material_variables[:,0] = default_n0   # n0
272material_variables[:,1] = default_d1   # d1 
273material_variables[:,2] = default_n1   # n1 
274material_variables[:,3] = default_d2   # d2 
275material_variables[:,4] = default_n2   # n2 
276
277print '         ++++ Initially - material_variables data set to default values across whole domain ' 
278
279# Update material_variables from model_data.mat_poly_list
280
281# Obtain a list of real world centroid coords of model elements
282points = domain.get_centroid_coordinates(absolute=True)   
283
284for k, poly in enumerate(model_data.mat_poly_list):
285   
286    # Get indices of triangles inside polygon k
287    indices = inside_polygon(points, poly, closed=True)
288
289    # Index into material data for polygon k
290    rid = model_data.mat_RID_list[k]
291       
292    # Material data for polygon     
293    n0, d1, n1, d2, n2 = model_data.mat_roughness_data[rid]   
294   
295    # Update the data for the elements inside the mat_poly   
296    for i in indices:
297        material_variables[i, 0] = n0
298        material_variables[i, 1] = d1       
299        material_variables[i, 2] = n1
300        material_variables[i, 3] = d2
301        material_variables[i, 4] = n2                       
302       
303   
304
305print '         ++++ Finally - material_variables data updated spatially with the material values set in the mif/mid files ' 
306
307   
308print model_data.basename+' >>>> Completed assignment of material (surface roughness) data to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
309
310print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
311
312
313
314#########################################################################################
315#                                                                                       #
316#                                     RUN THE SIMULATION                                #
317#                                                                                       #
318#########################################################################################
319
320#----------------------------------------------------------------------------------------
321# This section starts the analysis of flow through the domain from the specified (realworld)
322# starttime  up to the specified (realworld) finaltime
323# All time is in seconds.                                                 
324#
325# Note:
326# 1. -- The internal computational timestep is set at each timestep to meet a CFL=1 limit.
327#       It is therefore important to examine the mesh statistics to see that there are no
328#       unintended small triangles, particularly in areas of deeper water as they can
329#       dramatically reduce the computational timerstep and slow the model down.
330# 2. -- The model runs on 'internal' time which will normally differ from realworld or
331#       'Design' event time as input by the user in various datafiles. This 'internal'time is
332#       normally not seen by the user but users should be aware of its existence.
333# 3. -- In a similar fashion the model runs with all coordinates recomputed to a local
334#       lower left frame origin (to impove numerical precission). The output files are in these
335#       'internal' coordinates but each file also stores the new origin coordintes from which
336#       the realworld coordinates can be restored (as presently occurs in animate.exe)
337# 4. -- 'friction' is recomputed as a function of depth and location at each steptime for
338#       all wet cells. Recomputation at every computational step was considered unnecessay
339#       and greatly incresed run times.
340#----------------------------------------------------------------------------------------
341print model_data.basename+' >>>> Starting the simulation for ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event
342starttime     = 40*3600              # start simulation at t=40h relative to user/data time
343endtime       = starttime + 4*3600   # end simulation at t=44h relative to user/data time
344steptime      = 300                  # write sww  file results out every 300secs = 5min
345
346print model_data.basename+' >>>> Simulation is for %.2f hours with output to the sww file at %.0f minute intervals ' %(float((endtime-starttime)/3600.0),float(steptime/60.0))
347
348#  Initiate the starttime reset
349domain.set_starttime(starttime)
350print 'Start time', domain.get_time()
351
352
353for t in domain.evolve(yieldstep = steptime, finaltime = endtime):       # This is a steptime loop (not computaional timestep)
354   
355    domain.write_time()                                                  # confirm t each steptime
356    domain.write_boundary_statistics(tags='Lake', quantities='stage')    # indicate bdry stats each steptime
357    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
358
359    print '         ++++ Recomputing friction at t= ', t
360
361    # rebuild the 'friction' values adjusted for depth
362    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
363   
364        # Get roughness variables
365        n0 = material_variables[i,0]
366        d1 = material_variables[i,1]
367        n1 = material_variables[i,2]
368        d2 = material_variables[i,3]
369        n2 = material_variables[i,4]               
370               
371
372        # Recompute friction values from depth for this element
373        d = depth.get_values(location='centroids', indices=[i])[0]
374         
375        if d <= d1: 
376            depth_dependent_friction = n1
377        elif d >= d2:
378            depth_dependent_friction = n2                       
379        else:
380            depth_dependent_friction = n1+(n2-n1)/(d2-d1)*d                         
381
382
383        domain.set_quantity('friction',
384                            numeric=depth_dependent_friction,
385                            location='centroids',
386                            indices=[i],
387                            verbose=model_data.anuga_verbose)                               
388             
389             
390    if model_data.model_verbose :
391        friction = domain.get_quantity('friction')
392       
393        # Print something
394       
395       
396
397print model_data.basename+' >>>> ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event,' -- Simulation completed succesfully'
398print model_data.basename+' >>>> Completed ', endtime/3600.0, 'hours of model simulation at t= %.2f hours' %(float(time.time()-t0)/3600.0)
399
Note: See TracBrowser for help on using the repository browser.