Changeset 5989

Timestamp:

Nov 21, 2008, 3:35:51 PM (16 years ago)

Author:

ole

Message:

Update from Ted

File:

• anuga_work/development/ted_rigby/model_run.py (modified) (10 diffs)

anuga_work/development/ted_rigby/model_run.py

@@ -11 +11 @@
 
 Code Version : 1.00 June 2008 Initial release
+               1.01 December 2008 Vble_n recoded to speed execution
 Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au
 
@@ -141 +142 @@
 
 # Assign uniform initial friction value to domain mesh centroids as initial condition
-# Note need to re-assess initial friction value based on rough_polys read in from tuflow
-# files and computed depth during evolve time loop!!!!
+# Note: Friction will be re-assessed  based on roughness data read in from tuflow
+# 2d_mat files and computed depth during evolve time loop!!!!
 domain.set_quantity('friction', model_data.default_n, location = 'centroids')
 
@@ -201 +202 @@
 tms_filename=model_data.dsbc_hydrographs[0]
 lake_stage  = file_function(tms_filename,quantities='value',boundary_polygon=None)
-# Note: will need to be modified if more than one dsbcis to be applied
+# Note: will need to be modified if more than one dsbc is to be applied
 print model_data.basename + ' >>>> Applying dsbc %s to model at the Lake_bdry ' %(tms_filename)
 
@@ -227 +228 @@
 #######################################################################################
 #                                                                                     #
-#  READ IN THE (TEMPORALLY DEPTH DEPENDENT) SURFACE ROUGHNESS (FRICTION) PARAMETERS   #
+#  ASSIGN THE (TEMPORALLY DEPTH DEPENDENT) SURFACE ROUGHNESS (FRICTION) PARAMETERS    #
 #                                                                                     #
 #######################################################################################
 
 #---------------------------------------------------------------------------------------
-# This section creates a new domain quantity mat_RID and assigns surface roughness IDs (RID)
-# values to this quantity for later use in the evolve loop where RID is used with depth to
-# compute a depth dependent friction at each time step.
+# This section creates an array of surface roughness data for each element in the domain 
+# ready for later use in the evolve loop where the roughness data is  used with depth to
+# compute a depth dependent friction at each yieldstep (output timestep to the swww file).
 #
 # Note that because the polys read in should but may not cover the entire bounding poly area
-# It is necessarry to declare a default material as in Tuflow applying the default to mat_RID
-# across the whole bounding poly area before applying (overwriting) with the actual RID polys
-# and their RID values obtained from the 2d_mat mid/mif files!!!!
+# It is prudent to initially declare a default material (as in Tuflow) applying the default 
+# to the whole domain before applying (overwriting) with the spatial roughness data from the 
+# polys obtained from the 2d_mat mid/mif files!!!!
 #--------------------------------------------------------------------------------------
 
@@ -245 +246 @@
 
 N = len(domain)
-material_variables = zeros((N, 5), Float)
-
-# Set default material variables n0, d1, n1, d2, n2
-material_variables[:,0] = 0.035 # n0
-material_variables[:,1] = 0.3   # d1  
-material_variables[:,2] = 0.050 # n1  
-material_variables[:,3] = 1.5   # d2  
-material_variables[:,4] = 0.030 # n2  
-
-
-# Assign polygonal data to material_variables from model_data.mat_poly_list
-
-
-points = domain.get_centroid_coordinates(absolute=True)
+# Create the new roughness array to hold the elemental roughness data
+# In this array of Nx5 values
+#       n0 -- is the equivalent fixed roughness value  
+#       d1 -- is the depth below which n1 applies
+#       n1 -- is the roughness applying below d1
+#       d2 -- is the depth above which n2 applies
+#       n2 -- is the roughness applying above d2  
+
+# Create and assign zeros to the roughness array
+material_variables = zeros((N, 5), Float)                
+
+# Set default material variables n0, d1, n1, d2, n2 in the roughness array
+default_n0 = 0.040
+default_d1 = 0.3
+default_n1 = 0.050
+default_d2 = 1.5
+default_n2 = 0.030
+
+material_variables[:,0] = default_n0   # n0
+material_variables[:,1] = default_d1   # d1  
+material_variables[:,2] = default_n1   # n1  
+material_variables[:,3] = default_d2   # d2  
+material_variables[:,4] = default_n2   # n2  
+
+print '         ++++ Initially - material_variables data set to default values across whole domain ' 
+
+# Update material_variables from model_data.mat_poly_list
+
+# Obtain a list of real world centroid coords of model elements
+points = domain.get_centroid_coordinates(absolute=True)   
 
 for k, poly in enumerate(model_data.mat_poly_list):
@@ -270 +287 @@
     # Material data for polygon     
     n0, d1, n1, d2, n2 = model_data.mat_roughness_data[rid]    
-        
+    
+    # Update the data for the elements inside the mat_poly    
     for i in indices:
         material_variables[i, 0] = n0
@@ -280 +298 @@
     
 
-print '         ++++ Upgraded material_variables spatially with the material values set in the mif/mid files ' 
-
-    
-print model_data.basename+' >>>> Completed assignment of material (surface roughness) and default n to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
+print '         ++++ Finally - material_variables data updated spatially with the material values set in the mif/mid files ' 
+
+    
+print model_data.basename+' >>>> Completed assignment of material (surface roughness) data to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
 
 print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
@@ -296 +314 @@
 
 #----------------------------------------------------------------------------------------
-# This section starts the analysis of flow through the domain up to the specified finaltime
-# All time is in seconds                                                
+# This section starts the analysis of flow through the domain from the specified (realworld)
+# starttime  up to the specified (realworld) finaltime
+# All time is in seconds.                                                 
 #
-# Note: The internal computational timestep is set at each timestep to meet a CFL=1 limit.
-# It is therefore important to examine the mesh statistics to see that there are no 
-# unintended small triangles, particularly in areas of deeper water as DT=DX/(g*H)**0.5
-# Yieldstep is the interval at which results are saved to the output (sww) file
+# Note:
+# 1. -- The internal computational timestep is set at each timestep to meet a CFL=1 limit.
+#       It is therefore important to examine the mesh statistics to see that there are no 
+#       unintended small triangles, particularly in areas of deeper water as they can
+#       dramatically reduce the computational timerstep and slow the model down.
+# 2. -- The model runs on 'internal' time which will normally differ from realworld or
+#       'Design' event time as input by the user in various datafiles. This 'internal'time is
+#       normally not seen by the user but users should be aware of its existence.
+# 3. -- In a similar fashion the model runs with all coordinates recomputed to a local
+#       lower left frame origin (to impove numerical precission). The output files are in these
+#       'internal' coordinates but each file also stores the new origin coordintes from which
+#       the realworld coordinates can be restored (as presently occurs in animate.exe)
+# 4. -- 'friction' is recomputed as a function of depth and location at each steptime for
+#       all wet cells. Recomputation at every computational step was considered unnecessay
+#       and greatly incresed run times.
 #----------------------------------------------------------------------------------------
-
-# Note: 'friction' is recomputed as a function of depth and location at each timestep( now each steptime!!)
-
-# EHR - could reduce recalc of friction to say every steptime??? this would markedly speedup?? Done??
-# EHR - need to set a startime and use in domain.set_time(starttime) to restrict to area of data of interest
-#       not clear how initial conditions etc fit in to a start not at t=0???
-# EHR - could restrict update of friction to wet cells only as well
-
 print model_data.basename+' >>>> Starting the simulation for ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event
-starttime     = 40*3600            # start simulation at t=40h relative to data
-endtime       = 9000               # end at 9000 secs (150min) - Increase to 40hrs when debugged!!!!!!!!
-steptime      = 300                # write sww  file results out every 300secs = 5min
-
-print 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))
-
-#  
+starttime     = 40*3600              # start simulation at t=40h relative to user/data time
+endtime       = starttime + 4*3600   # end simulation at t=44h relative to user/data time
+steptime      = 300                  # write sww  file results out every 300secs = 5min
+
+print 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))
+
+#  Initiate the starttime reset
 domain.set_starttime(starttime)
 
 
-for t in domain.evolve(yieldstep = steptime, finaltime = endtime):
+for t in domain.evolve(yieldstep = steptime, finaltime = endtime):       # This is a steptime loop (not computaional timestep)
     
     domain.write_time()                                                  # confirm t each steptime
@@ -329 +351 @@
     depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
 
-
-    print '         ++++ Recomputing friction at t= ', domain.get_time()
-        
+    print '         ++++ Recomputing friction at t= ', t
 
     # rebuild the 'friction' values adjusted for depth 
-    for i in domain.get_wet_elements():                                        # loop for each wet element in domain
+    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
     
         # Get roughness variables 
@@ -344 +364 @@
                 
 
-        # Recompute friction values from depth for this element get_tuflow_data.py
+        # Recompute friction values from depth for this element 
         d = depth.get_values(location='centroids', indices=[i])[0]