Changeset 5139 for anuga_work/production/tonga/tongatapu.py

Timestamp:

Mar 7, 2008, 4:37:36 PM (15 years ago)

Author:

herve

Message:

deleting data files

File:

• anuga_work/production/tonga/tongatapu.py (modified) (3 diffs)

anuga_work/production/tonga/tongatapu.py

@@ -1 @@
-"""Script for running a tsunami inundation scenario for Cairns, QLD Australia.
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
 
 Source data such as elevation and boundary data is assumed to be available in
 directories specified by project.py
-The output sww file is stored in directory named after the scenario, i.e
-slide or fixed_wave.
+The output sww file is stored in project.output_run_time_dir 
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and a tsunami wave generated by a submarine mass failure.
-
-Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and
-Nick Bartzis, GA - 2006
+the elevation data and a simulated tsunami generated with URS code.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
 """
 
@@ -18 +16 @@
 
 # Standard modules
-import os
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
 import time
 import sys
@@ -24 +25 @@
 # Related major packages
 from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
 from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
-from anuga.shallow_water import Time_boundary
-from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Field_boundary
+from Numeric import allclose
+from anuga.shallow_water.data_manager import export_grid
+
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
-from anuga.shallow_water.data_manager import dem2pts
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
+from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+from anuga_parallel.parallel_abstraction import get_processor_name
+from anuga.caching import myhash
+from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
+from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
 
 # Application specific imports
 import project                 # Definition of file names and polygons
 
-
-#------------------------------------------------------------------------------
-# Define scenario as either slide or fixed_wave.
-#------------------------------------------------------------------------------
-scenario = 'slide' 
-#scenario = 'fixed_wave'
-
-if os.access(scenario, os.F_OK) == 0:
-    os.mkdir(scenario)
-basename = scenario + 'source'
-
-
-#------------------------------------------------------------------------------
-# Preparation of topographic data
-# Convert ASC 2 DEM 2 PTS using source data and store result in source data
-#------------------------------------------------------------------------------
-
-# Filenames
-dem_name = 'tongatapu_10mgrid' 
-meshname = 'tongatapu.msh'
-
-# Create DEM from asc data
-convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
-
-# Create pts file for onshore DEM
-dem2pts(dem_name, use_cache=True, verbose=True)
-
-
-#------------------------------------------------------------------------------
-# 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
-#------------------------------------------------------------------------------
-
-remainder_res = 1000000
-islands_res = 5000
-cairns_res = 5000 
-interior_regions = [[project.poly_tongatapu, tongatapu_res],
-                    [project.poly_island1, islands_res],
-                    [project.poly_island2, islands_res],
-                    [project.poly_island3, islands_res],
-                    [project.poly_island4, islands_res],
-                    [project.poly_island5, islands_res],
-                    [project.poly_island6, islands_res],
-                    [project.poly_island7, islands_res],
-                    [project.poly_island8, islands_res],
-                    [project.poly_island9, islands_res],
-                    [project.poly_island10, islands_res],
-                    [project.poly_island11, islands_res],
-                    [project.poly_island12, islands_res],
-                    [project.poly_island13, islands_res]]
-
-create_mesh_from_regions(project.bounding_polygon,
-                         boundary_tags={'top': [0],
-                                        'ocean_east': [1],
-                                        'bottom': [2],
-                                        'onshore': [3]},
-                         maximum_triangle_area=remainder_res,
-                         filename=meshname,
-                         interior_regions=interior_regions,
-                         use_cache=True,
-                         verbose=True)
-
-
-#------------------------------------------------------------------------------
-# 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(basename) 
-domain.set_datadir(scenario)
-domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-domain.set_minimum_storable_height(0.01)
-
-#print 'domain.tight_slope_limiters', domain.tight_slope_limiters
-domain.tight_slope_limiters = 0
-print 'domain.tight_slope_limiters', domain.tight_slope_limiters
-
-
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
-
-tide = 0.0
-domain.set_quantity('stage', tide)
-domain.set_quantity('friction', 0.0) 
-domain.set_quantity('elevation', 
-                    filename=dem_name + '.pts',
-                    use_cache=True,
-                    verbose=True,
-                    alpha=0.1)
-
-
-
-#------------------------------------------------------------------------------
-# Setup information for slide scenario (to be applied 1 min into simulation
-#------------------------------------------------------------------------------
-
-if scenario == 'slide':
-    # Function for submarine slide
-    from anuga.shallow_water.smf import slide_tsunami  
-    tsunami_source = slide_tsunami(length=35000.0,
-                                   depth=project.slide_depth,
-                                   slope=6.0,
-                                   thickness=500.0, 
-                                   x0=project.slide_origin[0], 
-                                   y0=project.slide_origin[1], 
-                                   alpha=0.0, 
-                                   domain=domain,
-                                   verbose=True)
-
-
-#------------------------------------------------------------------------------
-# Setup boundary conditions
-#------------------------------------------------------------------------------
-
-print 'Available boundary tags', domain.get_boundary_tags()
-
-Br = Reflective_boundary(domain)
-Bd = Dirichlet_boundary([tide,0,0])
-
-# 60 min square wave starting at 1 min, 50m high
-if scenario == 'fixed_wave':
+def run_model(**kwargs):
+    
+
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+    print "Processor Name:",get_processor_name()
+
+    #copy script must be before screen_catcher
+    #print kwargs
+
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+
+    # filenames
+#    meshes_dir_name = project.meshes_dir_name+'.msh'
+
+    # creates copy of code in output dir
+    print 'min triangles', project.trigs_min,
+    print 'Note: This is generally about 20% less than the final amount'
+
+    #--------------------------------------------------------------------------
+    # 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
+    #--------------------------------------------------------------------------
+
+    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
+    # causes problems with the ability to cache set quantity which takes alot of times
+    
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project.poly_all,
+                             boundary_tags=project.boundary_tags,
+                             maximum_triangle_area=project.res_poly_all,
+                             interior_regions=project.interior_regions,
+                             filename=project.meshes_dir_name+'.msh',
+                             use_cache=False,
+                             verbose=True)
+    barrier()
+
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #above don't work
+    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
+    print 'memory usage before del domain',mem_usage()
+       
+    print domain.statistics()
+    print 'triangles',len(domain)
+    
+    kwargs['act_num_trigs']=len(domain)
+    
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    if myid == 0:
+
+        print 'Setup initial conditions'
+
+        from polygon import Polygon_function
+        #following sets the stage/water to be offcoast only
+#        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
+#                                 geo_reference = domain.geo_reference)
+#        domain.set_quantity('stage', IC)
+        domain.set_quantity('stage',kwargs['tide'] )
+#        domain.set_quantity('stage', kwargs['tide'])
+        domain.set_quantity('friction', kwargs['friction']) 
+        
+        print 'Start Set quantity',kwargs['bathy_file']
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = True,
+                            verbose = True,
+                            alpha = kwargs['alpha'])
+        print 'Finished Set quantity'
+    barrier()
+    
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
+
+    if numprocs > 1:
+        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['aa_scenario_name'])
+    domain.set_datadir(kwargs['output_dir'])
+    domain.set_default_order(2) # Apply second order scheme
+    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
+    domain.set_store_vertices_uniquely(False)
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.tight_slope_limiters = 1
+    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+
+    #------------------------------------------------------------------------------
+    # Setup information for slide scenario (to be applied 1 min into simulation
+    #------------------------------------------------------------------------------
+
+    if scenario == 'slide':
+        # Function for submarine slide
+        from anuga.shallow_water.smf import slide_tsunami  
+        tsunami_source = slide_tsunami(length=35000.0,
+                                       depth=project.slide_depth,
+                                       slope=6.0,
+                                       thickness=500.0, 
+                                       x0=project.slide_origin[0], 
+                                       y0=project.slide_origin[1], 
+                                       alpha=0.0, 
+                                       domain=domain,
+                                       verbose=True)
+    
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
+
+    if scenario == 'fixed_wave':
     # FIXME (Ole): Change this to Transmissive Momentum as
     # suggested by Rajaraman (ticket:208)
-    Bw = Time_boundary(domain = domain,
-                       f=lambda t: [(60<t<3660)*50, 0, 0])
-    domain.set_boundary({'ocean_east': Bw,
-                         'bottom': Bd,
-                         'onshore': Bd,
-                         'top': Bd})
-
-# boundary conditions for slide scenario
-if scenario == 'slide':
-    domain.set_boundary({'ocean_east': Bd,
-                         'bottom': Bd,
-                         'onshore': Bd,
-                         'top': Bd})
-    
-
-#------------------------------------------------------------------------------
-# Evolve system through time
-#------------------------------------------------------------------------------
-
-import time
-t0 = time.time()
-
-from Numeric import allclose
-from anuga.abstract_2d_finite_volumes.quantity import Quantity
-
-if scenario == 'slide':
-    
-    for t in domain.evolve(yieldstep = 10, finaltime = 60): 
+        Bw = Time_boundary(domain = domain,
+                           f=lambda t: [(60<t<3660)*100, 0, 0])
+        domain.set_boundary({'ocean_east': Bw,
+                             'ocean_southeast': Bd,
+                             'land': Bd,
+                             'ocean_north': Bd,
+                             'ocean_west':})
+
+    # boundary conditions for slide scenario
+    if scenario == 'slide':
+        domain.set_boundary({'ocean_east': Bd,
+                             'ocean_southeast': Bd,
+                             'land': Bd,
+                             'ocean_north': Bd,
+                             'ocean_west':})
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']): 
         domain.write_time()
         domain.write_boundary_statistics(tags = 'ocean_east')     
-        
-    # add slide
-    thisstagestep = domain.get_quantity('stage')
-    if allclose(t, 60):
-        slide = Quantity(domain)
-        slide.set_values(tsunami_source)
-        domain.set_quantity('stage', slide + thisstagestep)
-
-    for t in domain.evolve(yieldstep = 10, finaltime = 5000, 
-                           skip_initial_step = True):
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean_east')
-
-if scenario == 'fixed_wave':
-
-    # save every two mins leading up to wave approaching land
-    for t in domain.evolve(yieldstep = 120, finaltime = 5000): 
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean_east')     
-
-    # save every 30 secs as wave starts inundating ashore
-    for t in domain.evolve(yieldstep = 10, finaltime = 10000, 
-                           skip_initial_step = True):
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean_east')
-            
-print 'That took %.2f seconds' %(time.time()-t0)
+
+    x, y = domain.get_maximum_inundation_location()
+    q = domain.get_maximum_inundation_elevation()
+
+    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+
+    print 'That took %.2f seconds' %(time.time()-t0)
+
+    #kwargs 'completed' must be added to write the final parameters to file
+    kwargs['completed']=str(time.time()-t0)
+    
+    if myid==0:
+        store_parameters(**kwargs)
+    barrier
+    
+    print 'memory usage before del domain1',mem_usage()
+
+def export_model(**kwargs):    
+    #store_parameters(**kwargs)
+    
+#    print 'memory usage before del domain',mem_usage()
+    #del domain
+    print 'memory usage after del domain',mem_usage()
+    
+    swwfile = kwargs['output_dir']+kwargs['scenario_name']
+    print'swwfile',swwfile
+    
+    export_grid(swwfile, extra_name_out = 'town',
+            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            timestep = None,
+            reduction = max,
+            cellsize = kwargs['export_cellsize'],
+            NODATA_value = -1E-030,
+            easting_min = project.eastingmin,
+            easting_max = project.eastingmax,
+            northing_min = project.northingmin,
+            northing_max = project.northingmax,
+            verbose = False,
+            origin = None,
+            datum = 'WGS84',
+            format = 'txt')
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    kwargs['est_num_trigs']=project.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project.host
+    kwargs['res_factor']=project.res_factor
+    kwargs['starttime']=project.starttime
+    kwargs['yieldstep']=project.yieldstep
+    kwargs['midtime']=project.midtime
+    kwargs['finaltime']=project.finaltime
+    
+    kwargs['output_dir']=project.output_run_time_dir
+    kwargs['bathy_file']=project.combined_dir_name+'.txt'
+#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
+    kwargs['file_name']=project.home+'detail.csv'
+    kwargs['aa_scenario_name']=project.scenario_name
+    kwargs['ab_time']=project.time
+    kwargs['res_factor']= project.res_factor
+    kwargs['tide']=project.tide
+    kwargs['user']=project.user
+    kwargs['alpha'] = project.alpha
+    kwargs['friction']=project.friction
+    kwargs['time_thinning'] = project.time_thinning
+    kwargs['dir_comment']=project.dir_comment
+    kwargs['export_cellsize']=project.export_cellsize
+    
+
+    run_model(**kwargs)
+      
+    if myid==0:
+        export_model(**kwargs)
+    barrier