Index: anuga_work/development/friction_UA_flume_2006/parallel_run_dam.py
===================================================================
--- anuga_work/development/friction_UA_flume_2006/parallel_run_dam.py	(revision 4048)
+++ anuga_work/development/friction_UA_flume_2006/parallel_run_dam.py	(revision 4048)
@@ -0,0 +1,185 @@
+"""Script for running a dam break simulation of U of Aberdeen's dam break tank.
+
+
+Ole Nielsen and Duncan Gray, GA - 2006 
+
+
+Issues
+If running this is hand-set-up parallel, the python files are overwritten.
+
+"""
+
+
+#----------------------------------------------------------------------------
+# Import necessary modules
+#----------------------------------------------------------------------------
+
+# Standard modules
+import time
+from time import localtime, strftime
+import sys
+from shutil import copy
+from os import path, sep
+from os.path import dirname  #, basename
+
+
+import pypar
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.abstract_2d_finite_volumes.region import Set_region 
+from anuga.fit_interpolate.interpolate import interpolate_sww2csv
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
+     copy_code_files
+
+
+# Scenario specific imports
+import project                 # Definition of file names and polygons
+import create_mesh
+
+
+def main(friction, outputdir_name=None, is_trial_run=False):
+    basename = 'zz' + str(friction)
+    if is_trial_run is True:
+        outputdir_name += '_test'
+        yieldstep = 1
+        finaltime = 31
+    else:
+        yieldstep = 0.01
+        finaltime = 31
+        
+
+    rank = pypar.rank()  #My id
+    if 0 == rank:
+        pro_instance = project.Project(['data','flumes','dam_2006'],
+                                        outputdir_name=outputdir_name)
+    else:
+        # wait for a bit..I don't know if I need this
+        # I do it to avoid crashes from the project code,
+        # which checks if I dir is created, and if it isn't
+        #creates it.
+        for i in range(500):
+            pass        
+        pro_instance = project.Project(['data','flumes','dam_2006'],
+                                       outputdir_name=outputdir_name)
+           
+    mesh_filename = pro_instance.meshdir + basename + '.msh'
+
+    #--------------------------------------------------------------------------
+    # Copy scripts to output directory and capture screen
+    # output to file
+    #--------------------------------------------------------------------------
+
+    # creates copy of code in output dir
+    print "The output dir is", pro_instance.outputdir
+    copy_code_files(pro_instance.outputdir,__file__,
+                    dirname(project.__file__) \
+                    + sep + project.__name__+'.py')
+    copy (pro_instance.codedir + 'run_dam.py',
+          pro_instance.outputdir + 'run_dam' + str(friction)+ '.py')
+    copy (pro_instance.codedir + 'create_mesh.py',
+          pro_instance.outputdir + 'create_mesh.py')
+    if is_trial_run is False:
+        start_screen_catcher(pro_instance.outputdir, int(friction*100))
+
+    print 'USER:    ', pro_instance.user
+    #-------------------------------------------------------------------------
+    # Create the triangular mesh
+    #-------------------------------------------------------------------------
+    
+    gate_position = 12.0
+    create_mesh.generate(mesh_filename,
+                         gate_position,
+                         is_coarse=is_trial_run) # this creates the mesh
+
+    head,tail = path.split(mesh_filename)
+    copy (mesh_filename,
+          pro_instance.outputdir + tail )
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    domain = Domain(mesh_filename, use_cache = False, verbose = True)
+   
+
+    print 'Number of triangles = ', len(domain)
+    print 'The extent is ', domain.get_extent()
+    print domain.statistics()
+
+    
+    domain.set_name(basename)
+    domain.set_datadir(pro_instance.outputdir)
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.set_minimum_storable_height(0.01)
+    #domain.set_store_vertices_uniquely(True)  # for writting to sww
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+
+    domain.set_quantity('stage', 0.0)
+    domain.set_quantity('friction', friction) 
+    domain.set_quantity('elevation',0.0)
+
+    domain.set_quantity('elevation', filename = 'elevation.xya',
+                        alpha = 10.0,verbose = True, use_cache = True)
+    
+    print 'Available boundary tags', domain.get_boundary_tags()
+    domain.set_region('dam','stage',0.6,
+                                 location = 'unique vertices') 
+    #domain.set_region(Set_region('dam','stage',0.2,
+    #                             location = 'unique vertices'))
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
+    domain.set_boundary( {'wall': Br, 'edge': Bd} )
+
+    #-------------------------------------------------------------------------
+    # Evolve system through time
+    #-------------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep, finaltime):
+        domain.write_time()
+    
+        print 'That took %.2f seconds' %(time.time()-t0)
+        print 'finished'
+
+    points = [[2.8,0.225],  #-1.8m from SWL
+              [5.1,0.225],  #0.5m from SWL
+              [6.6,0.225],  #2m from SWL
+              [6.95,0.255], #2.35m from SWL
+              [7.6,0.255],  #3m from SWL
+              [8.1,0.255],  #3.5m from SWL
+              [9.1,0.255]  #4.5m from SWL
+              ]
+
+#     points = [[gate_position - 0.65,0.2],
+#               [gate_position - 0.55,0.2],
+#               [gate_position - 0.45,0.2],
+#               [gate_position - 0.35,0.2],
+#               [gate_position - 0.25,0.2]
+#               ]
+
+    #-------------------------------------------------------------------------
+    # Calculate gauge info
+    #-------------------------------------------------------------------------
+
+    interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
+                        points,
+                        pro_instance.outputdir + "depth_manning_"+str(friction)+".csv",
+                        pro_instance.outputdir + "velocity_x.csv",
+                        pro_instance.outputdir + "velocity_y.csv")
+  
+    return pro_instance
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    rank = pypar.rank()  #My id
+    size = pypar.size() #Number of MPI processes
+    node = pypar.get_processor_name()
+    frictions = [0.0, 0.01, 0.0125]
+    for i, friction in enumerate(frictions):
+        if i%size == rank:
+            print "I am processor %d of %d on node %s" %(rank, size, node)
+            #main(friction, is_trial_run = True, outputdir_name='friction_UA')
