source: production/karratha_2005/run_karratha.py @ 1988

"""Script for running a tsunami inundation scenario for Karratha, 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 project.outputdir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and boundary data obtained from a tsunami simulation done with MOST.

Ole Nielsen, GA - 2005
"""

#tide = 0.75   #HMWS estimate by Colin French, GA
tide = 0       #HMW 


import os
import time


from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
     Dirichlet_boundary, Time_boundary, Transmissive_boundary
from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
     dem2pts, ferret2sww
from pyvolution.pmesh2domain import pmesh_to_domain_instance 
from caching import cache
import project

#Data preparation
#Convert ASC 2 DEM 2 PTS using source data and store result in source data
demname = project.demname

cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
      dependencies = [demname + '.asc'],
      verbose = True)
      #evaluate = True)

cache(dem2pts, demname, {'verbose': True},
      dependencies = [demname + '.dem'],      
      verbose = True)


#Convert MOST boundary
source_dir = project.boundarydir

from pyvolution.data_manager import ferret2sww

south = project.south 
north = project.north
west = project.west
east = project.east


cache(ferret2sww,
      (source_dir+project.boundary_basename,
       project.boundary_basename), 
      {'verbose': True,
       'minlat': south-1,
       'maxlat': north+1,
       'minlon': west-1,
       'maxlon': east+1,
       'origin': project.mesh_origin,
       'mean_stage': tide,
       'zscale': 1,                 #Enhance tsunami
       'fail_on_NaN': False,
       'inverted_bathymetry': True},
      #evaluate = True,
      verbose = True)
#FIXME: Dependencies

    
#ferret2sww(source_dir+project.boundary_basename,
#           project.boundary_basename, 
#           verbose=True,
#           minlat=south-1, maxlat=north+1,
#           minlon=west-1, maxlon=east+1,
#           origin = project.mesh_origin,
#           mean_stage = tide,
#           zscale = 1,
#           fail_on_NaN = False,
#           inverted_bathymetry = True)


#Create Triangular Mesh
from create_mesh import create_mesh

interior_regions = [[project.karratha_polygon, 25000],
                    [project.dampier_polygon, 8000],
                    [project.refinery_polygon, 8000]]

m = cache(create_mesh,
          project.polygon,
          {'boundary_tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]},
          'resolution': 100000,
          'filename': project.meshname + '.msh',
          'interior_regions': interior_regions},
          verbose = True)          
          #verbose = True,
          #evaluate = True )


#Setup domain
mesh = project.meshname + '.msh'
domain = cache(pmesh_to_domain_instance, (mesh, Domain),
               dependencies = [mesh],                     
               verbose = True)               


domain.set_name(project.basename)
domain.set_datadir(project.outputdir)
domain.store = True

domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
#domain.quantities_to_be_stored = ['stage']
        
print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()



#Setup Initial Conditions
domain.set_quantity('friction', 0)
domain.set_quantity('stage', tide)
domain.set_quantity('elevation',
                    filename = demname + '.pts',
                    use_cache = True,
                    verbose = True)



#Setup Boundary Conditions
print domain.get_boundary_tags()

Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
#domain.starttime = 3000  #Obtained from MOST
domain.starttime = 0  #Obtained from MOST

Br = Reflective_boundary(domain)
Bt = Transmissive_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])
Bw = Time_boundary(domain=domain,
                   f=lambda t: [(60<t<660)*4, 0, 0])


domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bf} ) #MOST tsunami

#domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bd} )  #Nothing


#Evolve
import time
t0 = time.time()

#for t in domain.evolve(yieldstep = 60, finaltime = 15000):
#    domain.write_time()
#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')

#for t in domain.evolve(yieldstep = 20, finaltime = 35000,
#                       skip_initial_step = True):
#    domain.write_time()
#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')

for t in domain.evolve(yieldstep = 60, finaltime = 40000):
                       #skip_initial_step = True):
    domain.write_time()
    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')       

print 'That took %.2f seconds' %(time.time()-t0)