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

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

Restricted main loop to wet elements

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