source: anuga_work/production/karratha_2005/run_karratha.py @ 5587

"""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 anuga.pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
     Dirichlet_boundary, Time_boundary, Transmissive_boundary
from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
     dem2pts, ferret2sww
from anuga.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 anuga.pyvolution.data_manager import ferret2sww

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


#Howard: Ignore this
#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 anuga.pmesh.create_mesh import create_mesh_from_regions

resolution = 4000
interior_regions = [#[project.karratha_polygon, 25000],
                    #[project.dampier_polygon, 2000],
                    #[project.refinery_polygon, 2000],
                    #[project.point_polygon, 2000]]
                    [project.neil1_polygon, resolution],
                    [project.neil2_polygon, 64000]]

#For visualisation
#interior_regions = [[project.wb_polygon, 400],
#                    [project.lwb_polygon, 8000]]

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


#Setup domain

#domain = cache(pmesh_to_domain_instance, (meshname, Domain),
#               dependencies = [meshname],                     
#               verbose = True)


#This is how it will be Howard!
domain = pmesh_to_domain_instance(meshname, Domain, use_cache=True)


domain.set_name(project.basename + '_%d' %resolution)
domain.set_datadir(project.outputdir)
domain.store = True
#domain.smooth = False


#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()


#File boundary.
#
#Takes a arbitrary simulation output in the form of an sww file and use it to provide a boundary.
#Values at every boundary point at every timestep are established by interpolation from the
#values in the sww boundary file.
#Boundary values at arbitrary timesteps needed in evolve are obtained by linear interpolation
#
#See also docstring for File_boundary in generic_boundary_conditions.py
#
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': Bd,'side': Bd, 'ocean': Bd} )  #Nothing


#Evolve
import time
t0 = time.time()

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

#for t in domain.evolve(yieldstep = 20, finaltime = 15000,
#                       skip_initial_step = True):
#    domain.write_time()
#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
    
#for t in domain.evolve(yieldstep = 10, finaltime = 40000,
#                       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)