Changeset 9354

Timestamp:

Oct 21, 2014, 3:30:43 PM (11 years ago)

Author:

davies

Message:

Updates to merwether validation test + a report

Location:

trunk/anuga_core

Files:

• source/anuga/utilities/plot_utils.py (modified) (2 diffs)
• validation_tests/case_studies/merewether/Observations (added)
• validation_tests/case_studies/merewether/Observations/ObservationPoints.csv (added)
• validation_tests/case_studies/merewether/Road (added)
• validation_tests/case_studies/merewether/Road/RoadPolygon.csv (added)
• validation_tests/case_studies/merewether/merewether.py (modified) (10 diffs)
• validation_tests/case_studies/merewether/plot_results.py (added)
• validation_tests/case_studies/merewether/produce_results.py (added)
• validation_tests/case_studies/merewether/project.py (modified) (1 diff)
• validation_tests/case_studies/merewether/results.tex (added)

trunk/anuga_core/source/anuga/utilities/plot_utils.py

@@ -1039 +1039 @@
 
 
-            if(bounding_polygon is not None):
+            if ( (bounding_polygon is not None) and (len(cut_points)>0)):
                 # Cut the points outside the bounding polygon
                 gridq[cut_points]= numpy.nan
@@ -1060 +1060 @@
     return
 
-def plot_triangles(p, adjustLowerLeft=False):
+def plot_triangles(p, adjustLowerLeft=False, values=None, color='k'):
     """ Add mesh triangles to a pyplot plot
-    """
+        
+       @param p = object holding sww vertex information (from util.get_output)
+       @param adjustLowerLeft = if TRUE, use spatial coordinates, otherwise use ANUGA internal coordinates     
+       @param values = list or array of length(p.x), or None. All triangles are assigned this value (for face plotting colors).
+       @param color = edge color for polygons (using matplotlib.colors notation)
+    """
+    import matplotlib
     from matplotlib import pyplot as pyplot
+    from matplotlib.collections import PolyCollection
     #
     x0=p.xllcorner
     x1=p.yllcorner 
-    #
+
+    # Get current figure information
+    #ax = pyplot.axes()
+
+    ##
+    vertices = []
     for i in range(len(p.vols)):
         k1=p.vols[i][0]
         k2=p.vols[i][1]
         k3=p.vols[i][2]
-        if(not adjustLowerLeft):
-            pyplot.plot([p.x[k1], p.x[k2], p.x[k3], p.x[k1]], [p.y[k1], p.y[k2], p.y[k3], p.y[k1]],'-',color='black')
+        if not adjustLowerLeft:
+            vertices.append([ [p.x[k1], p.y[k1]], [p.x[k2], p.y[k2]], [p.x[k3], p.y[k3]] ])
         else:
-            pyplot.plot([p.x[k1]+x0, p.x[k2]+x0, p.x[k3]+x0, p.x[k1]+x0], [p.y[k1]+x1, p.y[k2]+x1, p.y[k3]+x1, p.y[k1]+x1],'-',color='black')
-        #pyplot.plot([p.x[k3], p.x[k2]], [p.y[k3], p.y[k2]],'-',color='black')
-        #pyplot.plot([p.x[k3], p.x[k1]], [p.y[k3], p.y[k1]],'-',color='black')
+            vertices.append([ [p.x[k1]+x0, p.y[k1]+y0], [p.x[k2]+x0, p.y[k2]+y0], [p.x[k3]+x0, p.y[k3]+y0] ])
+      
+    if values is None: 
+        all_poly = PolyCollection( vertices, array = numpy.zeros(len(p.vols)), 
+            edgecolors=color)
+        all_poly.set_facecolor('none')
+    else:
+        assert len(values)==len(p.vols), 'len(values) must either be 1, or the same as len(p.vols)'
+        all_poly = PolyCollection( vertices, array = values, cmap = matplotlib.cm.jet, edgecolors=color)
+
+    pyplot.gca().add_collection(all_poly)
+    #    if(not adjustLowerLeft):
+    #        pyplot.plot([p.x[k1], p.x[k2], p.x[k3], p.x[k1]], [p.y[k1], p.y[k2], p.y[k3], p.y[k1]],'-',color=color)
+    #    else:
+    #        pyplot.plot([p.x[k1]+x0, p.x[k2]+x0, p.x[k3]+x0, p.x[k1]+x0], [p.y[k1]+x1, p.y[k2]+x1, p.y[k3]+x1, p.y[k1]+x1],'-',color=color)
+    #    #pyplot.plot([p.x[k3], p.x[k2]], [p.y[k3], p.y[k2]],'-',color='black')
+    #    #pyplot.plot([p.x[k3], p.x[k1]], [p.y[k3], p.y[k1]],'-',color='black')
 
 def find_neighbours(p,ind):

trunk/anuga_core/validation_tests/case_studies/merewether/merewether.py

@@ -7 +7 @@
 
 Water Research Laboratory, UNSW
+
+Edits by Steve Roberts & Gareth Davies, 2014
 """
 
@@ -22 +24 @@
 # Related major packages
 import anuga
-from anuga_parallel import myid, distribute, finalize
+from anuga_parallel import myid, distribute, finalize, barrier
 
 from anuga_parallel.parallel_operator_factory import Inlet_operator
+from anuga.utilities import quantity_setting_functions as qs
+from anuga import plot_utils as util
 
 # Application specific imports
 import project                 # Definition of file names and polygons
 
-verbose = project.verbose
+#verbose = project.verbose
 use_cache = project.use_cache
 
@@ -41 +45 @@
 meshname = 'merewether.msh'
 dem_name = 'topography1.dem'
+
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
 if myid == 0:
@@ -60 +68 @@
 holes_res = 1.0
 interior_regions = [[project.poly_merewether, merewether_res]]
-holes = project.holes
+
+# Either use houses as holes, or as part of the mesh (with breaklines)
+houses_as_holes = False
+if houses_as_holes:
+    holes = project.holes
+    breaklines = []
+else:
+    # Houses as elevation
+    house_height = 3.0
+    holes = []
+    breaklines = project.holes
+    house_addition_poly_fun_pairs = []
+    for i in range(len(breaklines)):
+        breaklines[i] = breaklines[i] + [breaklines[i][0]]
+        house_addition_poly_fun_pairs.append( 
+            [ breaklines[i], house_height])
+    house_addition_poly_fun_pairs.append(['All', 0.])
+
+#------------------------------------------------------------------------------
+# Make domain
+#------------------------------------------------------------------------------
 
 if myid == 0:
-    domain = anuga.create_domain_from_regions(project.bounding_polygon,
-                boundary_tags={'bottom': [0],
-                               'right': [1],
-                               'top': [2],
-                               'left': [3]},
-                                maximum_triangle_area=remainder_res,
-                                interior_holes=holes,
-                                mesh_filename=meshname,
-                                interior_regions=interior_regions,
-                                use_cache=use_cache,
-                                verbose=verbose)
+    mesh = anuga.create_mesh_from_regions(
+            project.bounding_polygon,
+            boundary_tags={'bottom': [0],
+                           'right': [1],
+                           'top': [2],
+                           'left': [3]},
+            maximum_triangle_area=remainder_res,
+            interior_holes=holes,
+            breaklines=breaklines,
+            filename=meshname,
+            interior_regions=interior_regions,
+            use_cache=use_cache,
+            verbose=verbose)
+
+    domain = anuga.create_domain_from_file(meshname)
+    
+    domain.set_flow_algorithm(alg)
 
     #------------------------------------------------------------------------------
@@ -79 +113 @@
     #------------------------------------------------------------------------------
     domain.set_quantity('stage', 0.0)
-    domain.set_quantity('friction', 0.02)
-    domain.set_quantity('elevation', filename='topography1.pts',
-                                          use_cache=use_cache,
-                                          verbose=verbose)
-
+
+    # Friction -- 2 options
+    variable_friction = True
+    if not variable_friction:
+        # Constant friction
+        domain.set_quantity('friction', 0.02)
+
+    else:
+        # Set friction to 0.02 on roads, 0.04 elsewhere
+        road_polygon = anuga.read_polygon('Road/RoadPolygon.csv')
+        friction_function = qs.composite_quantity_setting_function(
+            [ [road_polygon, 0.02], ['All', 0.04] ], 
+            domain)
+        domain.set_quantity('friction', friction_function)
+
+    # Elevation
+    if houses_as_holes:
+        domain.set_quantity('elevation', filename='topography1.pts',
+                              use_cache=use_cache,
+                                  verbose=verbose)
+
+    else:
+        domain.set_quantity('elevation', filename='topography1.pts',
+                              use_cache=use_cache,
+                              verbose=verbose, location='vertices')
+        # Add house_height inside houses    
+        house_addition_function = qs.composite_quantity_setting_function(
+            house_addition_poly_fun_pairs, domain)
+        domain.add_quantity('elevation', house_addition_function, 
+                            location='centroids')
 else:
     domain = None
@@ -93 +152 @@
 domain = distribute(domain)
 
+domain.set_store_vertices_uniquely()
+
 #------------------------------------------------------------------------------
 # Setup computational domain
@@ -98 +159 @@
 domain.set_name('merewether_1m') # Name of sww file
 domain.set_datadir('.') # Store sww output here
-domain.set_minimum_storable_height(0.001) # Store only depth > 1cm
+#domain.set_minimum_storable_height(0.001) # Store only depth > 1cm
 
 
@@ -112 +173 @@
 
 domain.set_boundary({'interior': Br,
-                     'bottom':   Br,
+                     'bottom':   Br,
                      'right':    Bt, # outflow
                      'top':      Bt, # outflow
@@ -121 +182 @@
 
 #fixed_inflow = anuga.Inflow(domain,
-#                       center=(382300.0,6354290.0),
-#                       radius=15.00,
-#                       rate=19.7)
+#           center=(382300.0,6354290.0),
+#           radius=15.00,
+#           rate=19.7)
 #domain.forcing_terms.append(fixed_inflow)
 #hydrograph = anuga.Inflow(center=(382300.0,6354290.0),radius=30.0,rate=anuga.file_function('test_hydrograph2.tms', quantities=['hydrograph']) 
@@ -132 +193 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep=10,finaltime=2000):
-        if myid == 0: print domain.timestepping_statistics()
+for t in domain.evolve(yieldstep=10,finaltime=1000):
+    if myid == 0: print domain.timestepping_statistics()
 
 
 domain.sww_merge()
 
+#------------------------------------------------------------------------------
+# Make geotif output
+#------------------------------------------------------------------------------
+
+barrier()
+if myid==0:
+    util.Make_Geotif('merewether_1m.sww', 
+                     output_quantities=['depth', 'velocity', 
+                                        'friction', 'elevation'],
+                     myTimeStep='last',
+                     CellSize=1.0,
+                     EPSG_CODE=32756,
+                     bounding_polygon=project.bounding_polygon,
+                     k_nearest_neighbours=1)
 
 finalize()

trunk/anuga_core/validation_tests/case_studies/merewether/project.py

@@ -19 +19 @@
 # bounding polygon for study area
 bounding_polygon = anuga.read_polygon('extent.csv')
-print 'Area of bounding polygon', anuga.polygon_area(bounding_polygon)/1000000.0
+print 'Area of bounding polygon (km^2)', anuga.polygon_area(bounding_polygon)/1.0e+06