Index: /anuga_work/development/gong_2008/project_slide.py
===================================================================
--- /anuga_work/development/gong_2008/project_slide.py	(revision 5675)
+++ /anuga_work/development/gong_2008/project_slide.py	(revision 5676)
@@ -20,4 +20,7 @@
 # INUNDATIONHOME is the inundation directory, not the data directory.
 home += sep +'data'
+
+# Bridgette to vary this parameter from 1 to 25 in steps of ?
+scale = None #25.
 
 #Making assumptions about the location of scenario data
@@ -53,5 +56,6 @@
 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
-outputtimedir = outputdir + local_time + sep
+dir_comment = '_scale_' + str(scale)
+outputtimedir = outputdir + local_time + dir_comment + sep
 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
 
Index: /anuga_work/development/gong_2008/project_slide_parameters.py
===================================================================
--- /anuga_work/development/gong_2008/project_slide_parameters.py	(revision 5676)
+++ /anuga_work/development/gong_2008/project_slide_parameters.py	(revision 5676)
@@ -0,0 +1,82 @@
+# -*- coding: cp1252 -*-
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+if sys.platform == 'win32':
+    home = getenv('INUNDATIONHOME')
+    user = getenv('USERPROFILE')
+
+else:    
+    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
+    user = getenv('LOGNAME')
+    print 'USER:', user
+
+###################################################################
+# Slide characteristics
+###################################################################
+
+# historical slides
+slide_origin_bulli = [370110,6189489]
+slide_origin_shovel = [354753,6187397]
+#slide_origin_yacaaba = [463012,6334735] original work was actually Birubi
+slide_origin_yacaaba = [476465, 6362424]
+slide_origin_birubi = [462267,6335100]
+
+# potential slides: these correspond approximately to modelled depths
+slide_origin_bulli_b = [405617,6249852]
+slide_origin_bulli_c0 = [438583,6308510]
+slide_origin_bulli_c1 = [465283,6335210]
+slide_origin_bulli_c2 = [483783,6358810]
+slide_origin_shovel_b = [396262,6252144]
+slide_origin_shovel_c0 = [431683,6311110]
+slide_origin_shovel_c1 = [458694,6340352]
+slide_origin_shovel_c2 = [475589,6363030]
+slide_origin_yacaaba_a = [359796,6190943]
+slide_origin_yacaaba_b = [397330,6251857]
+slide_origin_yacaaba_c0 = [431283,6312110]
+slide_origin_yacaaba_c1 = [459983,6339110]
+slide_origin_birubi_a = [362811,6191183]
+slide_origin_birubi_b = [403573,6250268]
+slide_origin_birubi_c0 = [437183,6309410]
+slide_origin_birubi_c2 = [483392,6359044]
+
+#bulli_depth = 2087.0
+bulli_depth = 1470.0
+bulli_length = 16840.0
+bulli_thickness = 424.0
+bulli_width = 8860.0
+bulli_density = 1.46
+bulli_slope = 4.0
+bulli_alpha = 126.0 - 90.0
+
+#shovel_depth = 968.0
+shovel_depth = 877.0
+shovel_length = 13500.0
+shovel_thickness = 165.0
+shovel_width = 4350.0
+shovel_density = 1.49
+shovel_slope = 4.0
+shovel_alpha = 118.0 - 90.0
+
+#yacaaba_depth = 1119.0
+yacaaba_depth = 938.
+yacaaba_length = 4189.
+yacaaba_thickness = 53.
+yacaaba_width = 2898.
+yacaaba_density = 1.48
+yacaaba_slope = 2.3
+yacaaba_alpha = 133.0 - 90.0
+ 
+birubi_depth = 1320.
+birubi_length = 9903.0
+birubi_width = 4150.
+birubi_density = 1.48
+birubi_slope = 3.7
+birubi_alpha = 133.0 - 90.0
+birubi_thickness = 140.
Index: /anuga_work/development/gong_2008/run_gong_slide_full.py
===================================================================
--- /anuga_work/development/gong_2008/run_gong_slide_full.py	(revision 5676)
+++ /anuga_work/development/gong_2008/run_gong_slide_full.py	(revision 5676)
@@ -0,0 +1,224 @@
+"""Script for running a tsunami inundation scenario for Wollongong, NSW, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_slide.py
+The output sww file is stored in project_slide.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by s submarine mass failure.
+
+Bridgette Lewis, 2008
+"""
+
+#-------------------------------------------------------------------------------
+# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os.path import dirname, basename
+from os import mkdir, access, F_OK, sep
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+
+# Application specific imports
+import project_slide              # Definition of file names and polygons
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+
+# creates copy of code in output dir 
+copy_code_files(project_slide.outputtimedir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
+myid = 0
+numprocs = 1
+start_screen_catcher(project_slide.outputtimedir, myid, numprocs)
+
+print 'USER:    ', project_slide.user
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#-------------------------------------------------------------------------------
+"""
+# filenames
+on_offshore10_dem_name = project_slide.on_offshore10_dem_name
+nsw_dem_name = project_slide.nsw_dem_name
+
+
+# creates DEM from asc data
+convert_dem_from_ascii2netcdf(on_offshore10_dem_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+dem2pts(on_offshore10_dem_name, use_cache=True, verbose=True)
+dem2pts(nsw_dem_name,
+        easting_min=project_slide.eastingmin_nsw,
+        easting_max=project_slide.eastingmax_nsw,
+        northing_min=project_slide.northingmin_nsw,
+        northing_max= project_slide.northingmax_nsw,
+        use_cache=True, verbose=True)
+
+print 'create offshore'
+G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.txt')
+G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.txt')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.txt')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.txt')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.txt')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.txt')+\
+      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.txt')
+print 'create onshore'
+G2 = Geospatial_data(file_name = project_slide.on_offshore10_dem_name + '.pts')
+G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
+print 'add'
+G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
+    G12.clip(Geospatial_data(project_slide.polyAll)) +\
+    G2.clip(Geospatial_data(project_slide.poly_10mclip)) +\
+    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_10mclip))
+print 'export points'
+G.export_points_file(project_slide.combined_dem_name + '.pts')
+#G.export_points_file(project_slide.combined_dem_name + '.xya')
+"""
+#----------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+meshname = project_slide.meshname+'.msh'
+remainder_res = 100000
+local_res = 25000
+gong_res = 500
+
+interior_regions = [[project_slide.poly_local, local_res],
+                    [project_slide.poly_gong, gong_res],
+                    [project_slide.poly_southgong, gong_res]]
+
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project_slide.polyAll,
+           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
+                              'e3': [3], 'e4': [4]},
+           'maximum_triangle_area': remainder_res,
+           'filename': meshname},
+           #'interior_regions': interior_regions},
+          verbose = True, evaluate=False)
+print 'created mesh'
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(meshname, use_cache = True, verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name(project_slide.basename)
+domain.set_datadir(project_slide.outputtimedir)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+#-------------------------------------------------------------------------------                                 
+# Set elevation to mesh
+#-------------------------------------------------------------------------------
+
+tide = 0.0
+domain.set_quantity('elevation',
+                    filename = project_slide.combined_dem_name + '.pts',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1
+                    )
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+from smf import slide_tsunami
+
+# effect on a3D and wavelength
+gamma = 1.85
+massco = 1.0
+dragco = 1.0
+
+# no effect on a3D and wavelength but used in Double Gaussian
+#dx = 0.01
+kappa = 3.
+kappad = 0.8
+
+# this doesn't seem to apper anywhere in smf
+frictionco = 0.01
+
+# Bridgette to vary this parameter from 1 to 25 in steps of (?) by changing it in the project file
+                  
+tsunami_source = slide_tsunami(project_slide.bulli_length,
+                               project_slide.bulli_depth,
+                               project_slide.bulli_slope,
+                               width=project_slide.bulli_width,
+			       thickness=project_slide.bulli_thickness, 
+                               x0=project_slide.slide_origin_bulli[0], 
+                               y0=project_slide.slide_origin_bulli[1], 
+                               alpha=project_slide.bulli_alpha,
+			       dx=None,			    
+ 			       scale=project_slide.scale,
+			       domain=domain,
+                               verbose=True)
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+domain.set_quantity('stage', tsunami_source)
+domain.set_quantity('friction', 0.0) 
+
+stage = domain.get_quantity('stage')
+etamax = stage.get_maximum_value()
+etamin = stage.get_minimum_value()
+
+filename = project_slide.outputtimedir + 'smf_data.csv'
+fid = open(filename,'w')
+s = 'scale, a3D, wavelength, etamax, etamin \n'
+fid.write(s)
+if project_slide.scale is not None:
+    s = '%.6f, %.6f, %.6f, %.6f, %.6f\n' \
+         %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
+else:
+    s = '%s, %.6f, %.6f, %.6f, %.6f\n' \
+        %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
+fid.write(s)   
+fid.close()
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+t0 = time.time()
+from Numeric import allclose
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
+    domain.write_time()
+        
+        
+print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'
Index: /anuga_work/development/gong_2008/run_gong_slide_parameters.py
===================================================================
--- /anuga_work/development/gong_2008/run_gong_slide_parameters.py	(revision 5676)
+++ /anuga_work/development/gong_2008/run_gong_slide_parameters.py	(revision 5676)
@@ -0,0 +1,96 @@
+"""Script for running a tsunami inundation scenario for Wollongong, NSW, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_slide_parameters.py
+The output sww file is stored in project_slide_parameters.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by s submarine mass failure.
+
+Bridgette Lewis, 2008
+"""
+
+#-------------------------------------------------------------------------------
+# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os.path import dirname, basename
+from os import mkdir, access, F_OK, sep
+import sys
+
+# Related major packages
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+from anuga.shallow_water import smf
+
+# Application specific imports
+import project_slide_parameters              # Definition of file names and polygons
+import project_slide
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+from smf import slide_tsunami
+
+# effect on a3D and wavelength
+length = project_slide_parameters.birubi_length
+width = project_slide_parameters.birubi_width
+depth = project_slide_parameters.birubi_depth
+slope = project_slide_parameters.birubi_slope
+thickness = project_slide_parameters.birubi_thickness
+alpha = project_slide_parameters.birubi_alpha
+x0 = project_slide_parameters.slide_origin_birubi[0]
+y0 = project_slide_parameters.slide_origin_birubi[1]
+
+gamma = 1.85
+massco = 1.0
+dragco = 1.0
+
+# no effect on a3D and wavelength but used in Double Gaussian
+dx = 0.01
+kappa = 3.
+kappad = 0.8
+
+# this doesn't seem to apper anywhere in smf
+frictionco = 0.01
+
+# scaling for Double Gaussian function
+scale = 30. # Bridgette's fiddle 
+
+fid = open('test.csv','w')
+s = 'parameter, wavelength, 3D_amp, eta_min\n'
+fid.write(s)
+
+for i in range(30):
+    scale = scale - 1.
+    #length = length + 10.
+    #depth = depth - 10.
+    tsunami_source = slide_tsunami(length=length,
+                                   width=width,
+                                   depth=depth,
+                                   slope=slope,
+                                   thickness=thickness, 
+                                   x0=x0, 
+                                   y0=y0, 
+                                   alpha=alpha,
+				   gamma=gamma,
+				   massco=massco,
+			           dragco=dragco,
+				   frictionco=frictionco,
+				   dx=dx,
+				   kappa=kappa,
+				   kappad=kappad,
+ 				   scale=scale,
+                                   verbose=False)
+
+    #parameter = depth
+    #parameter = length
+    parameter = scale
+
+    s = '%.6f, %.6f, %.6f, %.6f \n' %(parameter,tsunami_source.wavelength,tsunami_source.a3D,tsunami_source.nmin)
+    fid.write(s)
+
+fid.close()
