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

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

Removed duplicate lines in polys

File size: 22.2 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
13Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au
14
15-------------------------------------------------------------------------------
16"""
17
18#------------------------------------------------------------------------------
19# Import necessary modules
20#------------------------------------------------------------------------------
21
22# Standard python modules
23import os
24import time
25import sys
26
27# Related ANUGA modules
28from anuga.shallow_water import Domain
29from anuga.shallow_water import Reflective_boundary
30from anuga.shallow_water import Dirichlet_boundary
31from anuga.shallow_water import Time_boundary
32from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
33from anuga.shallow_water.shallow_water_domain import Inflow
34from anuga.abstract_2d_finite_volumes.util import file_function
35from anuga.abstract_2d_finite_volumes.quantity import Quantity
36from anuga.pmesh.mesh_interface import create_mesh_from_regions
37from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function
38from anuga.alpha_shape.alpha_shape import alpha_shape_via_files
39
40
41# Model specific imports
42import model_data                # Defines and provides the scenario-event specific data
43import get_tuflow_data           # Provides the routines to read in tuflow specific data
44
45#####################################################################################
46#                                                                                   #
47#                     S I M U L A T I O N    M A I N L I N E                        #
48#                                                                                   #
49#####################################################################################
50
51t0 = time.time()                 # record start time of this run (secs)
52print model_data.basename+' >>>> Comencing model_run.py code execution at t= %.2f hours'  %(float(time.time()-t0)/3600.0)
53
54#####################################################################################
55#                                                                                   #
56#                          CREATE THE TRIANGULAR MESH                               #
57#                                                                                   #
58#####################################################################################
59
60# Note: The mesh is created based on overall clipping polygon with a tagged 
61# boundary and interior regions as defined in model_data.py along with       
62# resolutions (maximal area per triangle) for each interior polygon       
63
64print model_data.basename+' >>>> Creating variable mesh within bounding_poly reflecting regional resolutions'
65# Create the mesh reflecting the bounding (default) and internal region mesh resolutions set in model_data.py
66# with boundary tags as defined in model_data.py
67
68create_mesh_from_regions(model_data.bounding_polygon,
69                         boundary_tags = model_data.boundary_strings,
70                         maximum_triangle_area = model_data.default_res,
71                         minimum_triangle_angle=30.0,
72                         interior_regions = model_data.interior_resregions,
73                         filename = model_data.mesh_filename,
74                         use_cache=True,
75                         verbose=model_data.anuga_verbose)
76
77
78
79print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
80
81######################################################################################
82#                                                                                    #
83#                       CREATE THE COMPUTATIONAL DOMAIN                              #
84#                                                                                    #
85######################################################################################
86
87print model_data.basename+' >>>> Creating domain from mesh and setting domain parameters'
88# create domain from mesh meshname
89domain = Domain(model_data.mesh_filename, use_cache=True, verbose=model_data.anuga_verbose) 
90print domain.statistics()                                # confirm what has been done
91
92# set base domain parameters
93domain.set_name(model_data.basename)                     # already created in model_data.py
94domain.set_datadir(model_data.results_dir)                 # already created and/or checked in model_data.py
95domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
96domain.set_minimum_storable_height(0.01)                 # Remove very shallow water depths from sww
97domain.set_maximum_allowed_speed(8.0)                    # Allow Maximum velocity of 8m/s....   
98
99######################################################################################
100#                                                                                    #
101#                       ASSIGN THE DOMAIN INITIAL CONDITIONS                         #
102#                                                                                    #
103######################################################################################
104
105#-------------------------------------------------------------------------------------
106# This section assigns the initial (t=0) conditions for the domain
107# Domain initial conditions can be fixed, read from a file, or from a function
108# Note also that anuga permits almost all initial conditions to subsequently vary
109# with time, so while the following conditions apply to the whole domain at t=0
110# all can be modified within the evolve time loop. In this way scour, deposition,
111# blockage, changes in roughness etc with time can be explicitly included in a simulation
112#-------------------------------------------------------------------------------------
113
114print model_data.basename+' >>>> Assigning initial conditions to %s domain' %(model_data.basename)
115
116# Elevation data may be assigned to the domain from a csv or (NetCDF) pts file.
117
118# For a smaller csv that will be read in repeatedly, best to use geospatial_data to readin and
119# export_point_file to save the points permanently in a NetCDF *.pts file. This runs faster
120# than the block read of a csv - but - is limited to about 6 million points.
121# Alternatively - by directly assigning elevation from a csv file can take advantage of the block
122# read in domain.set_quantity so the elevation dataset can be virtually unlimited in size.
123
124# In this simulation the size of the elevation dataset is to large for conversion to a (NetCDF)
125# pts format so data is read in from the raw csv file.
126
127domain.set_quantity('elevation', 
128                    filename = model_data.elev_filename,
129                    use_cache=True,
130                    verbose=model_data.anuga_verbose,
131                    alpha=0.01)
132
133
134# Set the intial water level (stage) across the domain             
135InitialWaterLevel = model_data.InitialWaterLevel             # At t=0 water surface is in this model at initial lake level (0.3m)
136if domain.get_quantity('elevation') < InitialWaterLevel:     # This can create ponds if land is below initial tailwater level
137    domain.set_quantity('stage', InitialWaterLevel)          # set initial stage at initial tailwater level (wet)
138else:       
139    domain.set_quantity('stage', expression='elevation')     # else set initial stage at land level (dry)
140
141
142# Assign uniform initial friction value to domain mesh centroids as initial condition
143# Note need to re-assess initial friction value based on rough_polys read in from tuflow
144# files and computed depth during evolve time loop!!!!
145domain.set_quantity('friction', model_data.default_n, location = 'centroids')
146
147
148
149
150#######################################################################################
151#                                                                                     #
152#             ASSIGN THE (TEMPORAL) SUBAREAL INFLOWS  TO THE DOMAIN                   #
153#                                                                                     #
154#######################################################################################
155
156#--------------------------------------------------------------------------------------
157# This section reads the local or total temporal inflows read in in model_data.py from
158# the Tuflow ts1 files and stored as NETCDF tms file in the tms_files subdirectory,
159# into the domain.
160#
161# Note: The current ANUGA inflow function 'pours' flow onto the domain within a circle
162# or polygonal area specified by the user at the specified temporal rate.
163# EHR -- the above inflows should eventually be augmented by  inflows at the boundary
164# (total as distinct from local hydrographs) to better correspond with conventional
165# flood modelling.
166#--------------------------------------------------------------------------------------
167
168print model_data.basename+" >>>> Assigning the inflow hydrographs from the (%s) tms files directory" %(model_data.tms_dir)
169
170
171inflow_fields=[]
172# read in each line of inflow_hydrographs[tms_filename,xcoord,ycoord] and assign as inflow hydrograph to domain
173for i in range(len(model_data.inflow_hydrographs)) :
174   
175    tms_filename = model_data.inflow_hydrographs[i][0]                   # extract the inflow hydro tms_filename
176    xcoord = model_data.inflow_hydrographs[i][1]                         # extract the inflow hydro xcoord
177    ycoord = model_data.inflow_hydrographs[i][2]                         # extract the inflow hydro ycoord
178    print model_data.basename+' >>>> Reading inflow hydrograph from %s at %7.2f %7.2f and appending to forcing terms ' %(tms_filename,xcoord,ycoord)
179       
180    # convert time series to temporal function
181    flowrate = file_function(tms_filename, quantities='value',boundary_polygon=None)
182    hydrograph = Inflow(domain,center=(xcoord,ycoord),radius=10, rate=flowrate )
183
184    # append hydrograph to domain
185    domain.forcing_terms.append(hydrograph)
186
187print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.inflow_hydrographs),float(time.time()-t0)/3600.0)
188
189
190#######################################################################################
191#                                                                                     #
192#        ASSIGN THE (TEMPORAL) DOWNSTREAM BOUNDARY CONDITIONS TO THE DOMAIN           #
193#                                                                                     #
194#######################################################################################         
195
196
197print model_data.basename+' >>>> Assigning the DSBC from the (%s) tms files directory ' %(model_data.tms_dir)
198print model_data.basename+' >>>> Available boundary tags are ', domain.get_boundary_tags()
199
200# convert dsbc time series to function of time
201tms_filename=model_data.dsbc_hydrographs[0]
202lake_stage  = file_function(tms_filename,quantities='value',boundary_polygon=None)
203# Note: will need to be modified if more than one dsbcis to be applied
204print model_data.basename + ' >>>> Applying dsbc %s to model at the Lake_bdry ' %(tms_filename)
205
206MaMt_bdry   = Reflective_boundary(domain)    # all reflective except lake as using inflow()
207Mriv_bdry   = Reflective_boundary(domain)    # Really redundant untill get boundary inflows
208Cbck_bdry   = Reflective_boundary(domain) 
209WFra_bdry   = Reflective_boundary(domain) 
210EFra_bdry   = Reflective_boundary(domain) 
211
212lake_bdry   = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=lake_stage)
213closed_bdry = Reflective_boundary(domain)                         
214
215# boundary conditions for current scenario
216domain.set_boundary({'Lake': lake_bdry,
217                     'MaMt': MaMt_bdry,
218                     'Mriv': Mriv_bdry,
219                     'Cbck': Cbck_bdry,
220                     'WFra': WFra_bdry,
221                     'EFra': EFra_bdry,
222                     'exterior': closed_bdry })  # this relates to the residue of the bounding poly not explicitly asigned above
223   
224print 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)
225
226
227#######################################################################################
228#                                                                                     #
229#  READ IN THE (TEMPORALLY DEPTH DEPENDENT) SURFACE ROUGHNESS (FRICTION) PARAMETERS   #
230#                                                                                     #
231#######################################################################################
232
233#---------------------------------------------------------------------------------------
234# This section creates a new domain quantity mat_RID and assigns surface roughness IDs (RID)
235# values to this quantity for later use in the evolve loop where RID is used with depth to
236# compute a depth dependent friction at each time step.
237#
238# Note that because the polys read in should but may not cover the entire bounding poly area
239# It is necessarry to declare a default material as in Tuflow applying the default to mat_RID
240# across the whole bounding poly area before applying (overwriting) with the actual RID polys
241# and their RID values obtained from the 2d_mat mid/mif files!!!!
242#--------------------------------------------------------------------------------------
243
244print model_data.basename+' >>>> Creating and assigning mat_RID (surface Roughness ID)instance for use in evolve'
245
246# Create an instance of a new domain quantity called mat_RID to store the surface Roughness IDs - fill initially with zeros
247mat_RID = Quantity(domain)
248if model_data.model_verbose :
249    print '         ++++ Created a zeroed domain quantity of  %s and length %i called mat_RID' %( type(mat_RID),len(mat_RID) )
250
251# Set  all of the RID values initially to one fixed default value (makes sure all elements have a RID)
252mat_RID.set_values(numeric=model_data.default_RID,location='centroids',verbose=True)
253if model_data.model_verbose :
254    print '         ++++ Initially filled mat_RID with the specified default RID value of %i ' %( model_data.default_RID )
255
256# Assign  model_data.mat_RID_list values to domain.from model_data.mat_poly_list
257for i in range(len(model_data.mat_RID_list)) :
258    mat_RID.set_values(numeric=model_data.mat_RID_list[i],polygon=model_data.mat_poly_list[i])
259if model_data.model_verbose :
260    print '         ++++ Upgraded mat_RID spatially with the values set in the mif/mid files ' 
261
262
263# Dump some RIDS to confirm assignment occured as intended
264if model_data.model_verbose :
265    print '         ++++ Check printout of RID assignments to centroids of elements nos 100 200 and 300'
266    print '         ++++ Element 100s RIDs ', mat_RID.get_values(location='centroids',indices=[100])
267    print '         ++++ Element 200s RIDs ', mat_RID.get_values(location='centroids',indices=[200])
268    print '         ++++ Element 300s RIDs ', mat_RID.get_values(location='centroids',indices=[300])
269   
270print model_data.basename+' >>>> Completed assignment of material (surface roughness) RIDs and default n to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
271
272print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
273
274
275#########################################################################################
276#                                                                                       #
277#                                     RUN THE SIMULATION                                #
278#                                                                                       #
279#########################################################################################
280
281#----------------------------------------------------------------------------------------
282# This section starts the analysis of flow through the domain up to the specified finaltime
283# All time is in seconds                                               
284#
285# Note: The internal computational timestep is set at each timestep to meet a CFL=1 limit.
286# It is therefore important to examine the mesh statistics to see that there are no
287# unintended small triangles, particularly in areas of deeper water as DT=DX/(g*H)**0.5
288# Yieldstep is the interval at which results are saved to the output (sww) file
289#----------------------------------------------------------------------------------------
290
291# Note: 'friction' is recomputed as a function of depth and location at each timestep( now each steptime!!)
292
293# EHR - could reduce recalc of friction to say every steptime??? this would markedly speedup?? Done??
294# EHR - need to set a startime and use in domain.set_time(starttime) to restrict to area of data of interest
295#       not clear how initial conditions etc fit in to a start not at t=0???
296# EHR - could restrict update of friction to wet cells only as well
297
298print model_data.basename+' >>>> Starting the simulation for ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event
299starttime     = 0                  # start simulation at t=0 EHR not yet implemented ???????
300endtime       = 9000               # end at 9000 secs (150min) - Increase to 40hrs when debugged!!!!!!!!
301steptime      = 300                # write sww  file results out every 300secs = 5min
302lastprintstep = 0                  # initialise lastprintstep
303print model_data.basename+' >>>> Simulation is for %.2f hours with output to the sww file at %.0f minute intervals ' %(float(endtime/3600.0),float(steptime/60.0))
304
305
306for t in domain.evolve(yieldstep = steptime, finaltime = endtime):
307   
308    domain.write_time()                                                            # confirm t each steptime
309    domain.write_boundary_statistics(tags='Lake', quantities='stage')              # indicate bdry stats each steptime
310    depth = domain.create_quantity_from_expression('stage - elevation')            # create depth instance for this timestep
311
312    if int(float((t+0.00001)/steptime)) > lastprintstep  :                                   # Recompute friction at each print timestep
313        lastprintstep = int(float((t+0.00001)/steptime))
314        print '         ++++ Recomputing friction at t= ', t
315       
316        noelementsupdated = 0                                                      # intialise update check
317        # rebuild the 'friction' values adjusted for depth
318        for i in range(len(domain)) :                                              # loop for each element in domain
319       
320            element_RID_value=mat_RID.get_values(location='centroids',indices=[i]) # get current elements centroidal RID value
321        #    print '         ++++ Element %i centroidal element_RID_value is %s' %( i,element_RID_value[0] )
322            for k in range(len(model_data.mat_roughness_data_list)):               # loop for each matl in tmf list
323                if int(element_RID_value[0]) == int( model_data.mat_roughness_data_list[k][0] )  :  # got tmf line for that matl
324        #            print '         ++++ found tmf match for RID %i of %s ' %( int(element_RID_value[0]),model_data.mat_roughness_data_list[k][8].rstrip(' ')  )
325                    d1 = float(model_data.mat_roughness_data_list[k][4])
326                    n1 = float(model_data.mat_roughness_data_list[k][5])
327                    d2 = float(model_data.mat_roughness_data_list[k][6])
328                    n2 = float(model_data.mat_roughness_data_list[k][7])
329                    noelementsupdated = noelementsupdated +1                        # increment update counter start at zeroth element
330                    if   depth.get_values(location='centroids',indices=[i]) <= d1 : # recompute friction values from depth for this element get_tuflow_data.py
331
332                        domain.set_quantity('friction',numeric=n1,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
333                    elif depth.get_values(location='centroids',indices=[i]) >= d2 :
334                        domain.set_quantity('friction',numeric=n2,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
335                    else :
336                        n3 = n1+(n2-n1)/(d2-d1)* ( depth.get_values(location='centroids',indices=[i])[0] )
337                        domain.set_quantity('friction',numeric=n3,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
338
339        if i <> (noelementsupdated-1) :                                              # Check all elements actually updated
340            # EHR - will need to modify if only wet elements updated
341            print '         ++++ WARNING - %i elements in domain but only %i updated ' %(i+1,noelementsupdated)
342             
343    if model_data.model_verbose :
344        print "         ++++  Check  some  RID,depth, friction values returned during this timestep "   
345        print "               element 100 ",mat_RID.get_values(location='centroids',indices=[100]),depth.get_values(location='centroids',indices=[100]),domain.get_quantity('friction',location='centroids',indices=[100] )
346        print "               element 200 ",mat_RID.get_values(location='centroids',indices=[200]),depth.get_values(location='centroids',indices=[200]),domain.get_quantity('friction',location='centroids',indices=[200] )
347        print "               element 300 ",mat_RID.get_values(location='centroids',indices=[300]),depth.get_values(location='centroids',indices=[300]),domain.get_quantity('friction',location='centroids',indices=[300] )
348                       
349 
350
351       
352
353print model_data.basename+' >>>> ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event,' -- Simulation completed succesfully'
354print model_data.basename+' >>>> Completed ', endtime/3600.0, 'hours of model simulation at t= %.2f hours' %(float(time.time()-t0)/3600.0)
355
Note: See TracBrowser for help on using the repository browser.