Changeset 3883

Timestamp:

Oct 30, 2006, 11:41:12 AM (19 years ago)

Author:

ole

Message:

More cleanup of okushiri prior to release

Location:

anuga_validation/okushiri_2005

Files:

• create_okushiri.py (modified) (9 diffs)
• extract_timeseries.py (modified) (1 diff)
• project.py (modified) (3 diffs)
• run_okushiri.py (modified) (2 diffs)

anuga_validation/okushiri_2005/create_okushiri.py

@@ -8 +8 @@
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.geospatial_data import Geospatial_data
 
 import project
 
+def prepare_bathymetry(filename):
+    """Convert benchmark 2 bathymetry to NetCDF pts file.
+    This is a 'throw-away' code taylor made for files like
+    'Benchmark_2_bathymetry.txt' from the LWRU2 benchmark
+    """
+
+    print 'Creating', filename
+    
+    # Read the ascii (.txt) version of this file,
+    # make it comma separated and invert the bathymetry
+    # (Below mean sea level should be negative)
+    infile = open(filename[:-4] + '.txt')
+
+    points = []
+    attribute = []
+    for line in infile.readlines()[1:]: #Skip first line (the header)
+        fields = line.strip().split()
+
+        x = float(fields[0])
+        y = float(fields[1])
+        z = -float(fields[2]) # Bathymetry is inverted in original file
+        
+        points.append([x,y])
+        attribute.append(z)
+    infile.close()
+
+    # Convert to geospatial data and store as NetCDF
+    G = Geospatial_data(data_points=points,
+                        attributes=attribute)
+    G.export_points_file(filename)
+    
+
 
 def prepare_timeboundary(filename):
-    """Converting benchmark 2 time series to NetCDF tms file.
+    """Convert benchmark 2 time series to NetCDF tms file.
     This is a 'throw-away' code taylor made for files like
     'Benchmark_2_input.txt' from the LWRU2 benchmark
     """
 
-    textversion = filename[:-4] + '.txt'
-
-    print 'Preparing time boundary from %s' %textversion
+    from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import array
 
-    fid = open(textversion)
 
-    #Skip first line
+    print 'Creating', filename
+
+    # Read the ascii (.txt) version of this file
+    fid = open(filename[:-4] + '.txt')
+
+    # Skip first line
     line = fid.readline()
 
-    #Read remaining lines
+    # Read remaining lines
     lines = fid.readlines()
     fid.close()
@@ -44 +79 @@
 
 
-    #Create tms file
-    from Scientific.IO.NetCDF import NetCDFFile
+    # Create tms NetCDF file
 
-    print 'Writing to', filename
     fid = NetCDFFile(filename, 'w')
-
     fid.institution = 'Geoscience Australia'
     fid.description = 'Input wave for Benchmark 2'
@@ -73 +105 @@
 
 
-    #Preparing points
-    #(Only when using original Benchmark_2_Bathymetry.txt file)
-    # FIXME: Replace by using geospatial data
-    #
-    #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
-    #xya2pts(project.bathymetry_filename, verbose = True,
-    #        z_func = lambda z: -z)
-
+    # Prepare bathymetry
+    prepare_bathymetry(project.bathymetry_filename)
 
     # Prepare time boundary
@@ -94 +120 @@
 
 
-    base_resolution = 1
+    base_resolution = 1 # Use this to coarsen or refine entire mesh. 
 
-    # Basic geometry
+    # Basic geometry and bounding polygon
     xleft   = 0
     xright  = 5.448
@@ -102 +128 @@
     ytop    = 3.402
 
-    # Outline 
     point_sw = [xleft, ybottom]
     point_se = [xright, ybottom]
@@ -108 +133 @@
     point_ne = [xright, ytop]
 
-    # Midway points (left)
-    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
-    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
+    bounding_polygon = [point_se,
+                        point_ne,
+                        point_nw,
+                        point_sw]
+    
 
-    #Localised refined area for gulleys
+    # Localised refined area for gulleys
     xl = 4.8
     xr = 5.3
@@ -121 +148 @@
     p2 = [xr, yt]
     p3 = [xl, yt]
+    
     gulleys = [p0, p1, p2, p3]
+    
 
-    #Island area
+    # Island area and drawdown region
     island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
     island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
@@ -132 +161 @@
     island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
     island_7 = [xl-.01, yt]
-
-    
-
+ 
     island = [island_0, island_1, island_2,
               island_3, island_4, island_5,
@@ -140 +167 @@
 
 
-    #Midway points just inside boundary and associated rectangular region
-    point_mtop0 = [xleft + (xright-xleft)/3+1, ytop-0.02]
-    point_mbottom0 = [xleft + (xright-xleft)/3+1, ybottom+0.02]
-    point_se0 = [xright-0.02, ybottom+0.02]
-    point_ne0 = [xright-0.02, ytop-0.02]    
+    # Region spanning half right hand side of domain just inside boundary
+    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
+    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
+    rhs_se = [xright-0.02, ybottom+0.02]
+    rhs_ne = [xright-0.02, ytop-0.02]    
 
-    mid = [point_mtop0, point_ne0, point_se0, point_mbottom0]
+    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
 
 
-    # Bounding polygon and interior regions
-    bounding_polygon = [point_se,
-                        point_ne,
-                        point_mtop,
-                        point_nw,
-                        point_sw,
-                        point_mbottom]
     
-
-    interior_regions = [[mid, 0.0005],
+    # Interior regions and creation of mesh
+    interior_regions = [[rhs_region, 0.0005],
                         [gulleys, 0.00002*base_resolution],
                         [island, 0.00007*base_resolution]] 
 
-    
-
-    print 'start create mesh from regions'
     meshname = project.mesh_filename + '.msh'
     m = create_mesh_from_regions(bounding_polygon,
-                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
-                                                'wave': [3]},
+                                 boundary_tags={'wall': [0, 1, 3],
+                                                'wave': [2]},     
                                  maximum_triangle_area=0.1*base_resolution,
                                  interior_regions=interior_regions,

anuga_validation/okushiri_2005/extract_timeseries.py

@@ -56 +56 @@
 
 
-expected_covariances = {'Boundary': 5.288392008865989237e-05, 
-                        'ch5': 1.166748190444680592e-04,
-                        'ch7': 1.121816242516757850e-04,
-                        'ch9': 1.249543278366777640e-04}
+#expected_covariances = {'Boundary': 5.288392008865989237e-05, 
+#                        'ch5': 1.166748190444680592e-04,
+#                        'ch7': 1.121816242516757850e-04,
+#                        'ch9': 1.249543278366777640e-04}
+#
+#expected_differences = {'Boundary':  8.373150808730501615e-04,
+#                        'ch5': 3.425914311580337875e-03,
+#                        'ch7': 2.802327594773105189e-03,
+#                        'ch9': 3.198733498646373370e-03}
 
-expected_differences = {'Boundary':  8.373150808730501615e-04,
-                        'ch5': 3.425914311580337875e-03,
-                        'ch7': 2.802327594773105189e-03,
-                        'ch9': 3.198733498646373370e-03}
+
+expected_covariances = {'Boundary': 5.278365381306500051e-05, 
+                        'ch5': 1.168250251701857399e-04,
+                        'ch7': 1.125146958894896555e-04,
+                        'ch9': 1.248634476898675049e-04}
+
+expected_differences = {'Boundary': 8.502104901511024380e-04, 
+                        'ch5': 3.442302747303998579e-03,
+                        'ch7': 2.802987591027187534e-03,
+                        'ch9': 3.191368834008508938e-03}
+
+
 
 

anuga_validation/okushiri_2005/project.py

@@ -1 +1 @@
 """Common filenames for Okushiri Island validation
+Formats are given as ANUGA native netCDF where applicable.
+
 """
 
@@ -8 +10 @@
 validation_filename = 'output_ch5-7-9.txt'
 
-# A simple conversion from txt file to the more compact NetCDF format
+# Digital Elevation Model
 bathymetry_filename = 'Benchmark_2_Bathymetry.pts'
 
@@ -15 +17 @@
 
 # Model output
-#output_filename = 'okushiri_old_limiters.sww'
-#output_filename = 'okushiri_new_limiters.sww'
-#output_filename = 'okushiri_nl_coarse10_mesh.sww'
-#output_filename = 'okushiri_nl_coarse100_mesh.sww'
-#output_filename = 'okushiri_original_mesh.sww'
-#utput_filename = 'okushiri_new_mesh.sww'
-output_filename = 'okushiri_new_mesh1.sww'
+output_filename = 'okushiri.sww'
 
+
+

anuga_validation/okushiri_2005/run_okushiri.py

@@ -38 +38 @@
 print domain.statistics()
 
+
 #-------------------------
 # Initial Conditions
@@ -49 +50 @@
                     use_cache=True)
 
+
 #-------------------------
-# Distribute domain
+# Distribute domain if run in parallel
 #-------------------------
-domain = distribute(domain)
+if numprocs > 1:
+    domain = distribute(domain)
 
-# Parameters
+
+#-------------------------
+# Set simulation parameters
+#-------------------------
 domain.set_name(project.output_filename)  # Name of output sww file
 domain.set_default_order(2)               # Apply second order scheme