Changeset 3601 for anuga_work/production/MOST_example/MOST_example

Timestamp:

Sep 15, 2006, 12:30:30 PM (18 years ago)

Author:

duncan

Message:

 

Location:

anuga_work/production/MOST_example/MOST_example

Files:

• project.py (modified) (4 diffs)
• run_pt_hedland.py (modified) (9 diffs)

anuga_work/production/MOST_example/MOST_example/project.py

@@ -7 +7 @@
 #from anuga.utilities.polygon import read_polygon
 import sys
-from anuga.pmesh.create_mesh import convert_points_from_latlon_to_utm
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
 from time import localtime, strftime
@@ -52 +51 @@
 #Derive subdirectories and filenames
 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
-outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep
+outputtimedir = time+sep
 
 #print 'outputtimedir', outputtimedir
@@ -80 +79 @@
 tidal_filename = tidedir + 'pt_hedland_tide.txt'
 
-meshname = meshdir + basename
+meshname = basename+time
+mesh_elevname = basename+time+'_elev'
 #onshore_dem_name = datadir + onshore_name_dted
 onshore_dem_name = datadir + onshore_name_dli
@@ -86 +86 @@
 offshore_dem_name2 = datadir + offshore_name2
 coast_dem_name = datadir + coast_name
-combined_dem_name = datadir + 'pt_hedland_combined_elevation'
+combined_dem_name = 'pt_hedland_combined_elevation'
 outputname = outputtimedir + basename  #Used by post processing
 

anuga_work/production/MOST_example/MOST_example/run_pt_hedland.py

@@ -31 +31 @@
 from anuga.pyvolution.util import Screen_Catcher
 
+from anuga.fit_interpolate.fit import fit_to_mesh_file
+
 # Application specific imports
 import project                 # Definition of file names and polygons
@@ -47 +49 @@
 
 # normal screen output is stored in 
-screen_output_name = project.outputtimedir + "screen_output.txt"
-screen_error_name = project.outputtimedir + "screen_error.txt"
+#screen_output_name = project.outputtimedir + "screen_output.txt"
+#screen_error_name = project.outputtimedir + "screen_error.txt"
 
 # used to catch screen output to file
-sys.stdout = Screen_Catcher(screen_output_name)
-sys.stderr = Screen_Catcher(screen_error_name)
+#sys.stdout = Screen_Catcher(screen_output_name)
+#sys.stderr = Screen_Catcher(screen_error_name)
 print 'USER:    ', project.user
 
@@ -64 +66 @@
 
 # filenames
-meshname = project.meshname+'.msh'
+meshname = project.meshname + '.tsh'
+mesh_elevname = project.mesh_elevname  + '.tsh'
 source_dir = project.boundarydir
 
-# fine data (clipping the points file to smaller area)
-# creates DEM from asc data 
-convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True)
-
-#creates pts file from DEM
-dem2pts(project.onshore_dem_name,
-        easting_min=project.eastingmin,
-        easting_max=project.eastingmax,
-        northing_min=project.northingmin,
-        northing_max= project.northingmax,
-        use_cache=True,
-        verbose=True)
-
-print 'create G1'
-G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
-print 'create G2'
-G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
-print 'create G3'
-G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
-print 'create G4'
-G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
-print 'add G1+G2+G3+G4'
-G = G1 + G2 + G3 + G4
-print 'export G'
-G.export_points_file(project.combined_dem_name + '.pts')
 
 #-------------------------------------------------------------------------------                                 
@@ -101 +79 @@
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 
-region_res = 500000
-coast_res = 500
-pt_hedland_res = 5000
+region_res = 5000000
+coast_res = 500000
+pt_hedland_res = 500000
 interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
                     [project.poly_region, region_res]]
 
 print 'number of interior regions', len(interior_regions)
-
-from anuga.utilities.polygon import plot_polygons
-if sys.platform == 'win32':
-    #figname = project.outputtimedir + 'pt_hedland_polys'
-    figname = 'pt_hedland_polys_test'
-    plot_polygons([project.polyAll,project.poly_pt_hedland,project.poly_region],
-              figname,
-              verbose = True)    
 
 print 'start create mesh from regions'
@@ -125 +95 @@
                              'bottom1': [4], 'bottom2': [5],
                              'bottom3': [6], 'right': [7]},
-           'maximum_triangle_area': 250000,
+           'maximum_triangle_area': 5000000,
            'filename': meshname,           
            'interior_regions': interior_regions},
           verbose = True, evaluate=True)
 
+cache(fit_to_mesh_file,(meshname,
+                 'pt_hedland_combined_elevation_31204' + '.pts',
+                 mesh_elevname),
+      {'verbose': True}
+      #,evaluate = True      
+      ,verbose = False
+      )
 #-------------------------------------------------------------------------------                                 
 # Setup computational domain
 #-------------------------------------------------------------------------------                                 
-domain = Domain(meshname, use_cache = False, verbose = True)
+domain = Domain(mesh_elevname, use_cache = False, verbose = True)
 
 print domain.statistics()
@@ -158 +135 @@
 print 'hi and file',project.combined_dem_name + '.pts'
 
-domain.set_quantity('elevation', 
-                    filename = project.combined_dem_name + '.pts',
-                    use_cache = True,
-                    verbose = True,
-                    alpha = 0.1
-                    )
+#domain.set_quantity('elevation', 
+#                    filename = project.combined_dem_name + '.pts',
+#                    use_cache = True,
+#                    verbose = True,
+#                    alpha = 0.1
+ #                   )
 
 #-------------------------------------------------------------------------------                                 
@@ -180 +157 @@
 #note only need to do when an SWW file for the MOST boundary doesn't exist
 cache(ferret2sww,
-      (source_dir + project.boundary_basename,
-       source_dir + project.boundary_basename+'_'+project.basename), 
+      (project.boundary_basename,
+       project.boundary_basename+'_'+project.basename), 
       {'verbose': True,
        'minlat': south,
@@ -190 +167 @@
        'origin': domain.geo_reference.get_origin(),
        'mean_stage': tide,
-       'zscale': 1,                 #Enhance tsunami
+       'zscale': 10,                 #Enhance tsunami
        'fail_on_NaN': False,
        'inverted_bathymetry': True},
@@ -199 +176 @@
 print 'Available boundary tags', domain.get_boundary_tags()
 
-Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
+Bf = File_boundary(project.boundary_basename+'_'+project.basename + '.sww', 
                     domain, verbose = True)
 Br = Reflective_boundary(domain)