source: production/gippsland_2005/run_small_gippsland.py @ 2058

"""
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.

Duncan Gray, Ole Neilsen, GA - 2005
"""


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, asc_csiro2sww
#from pyvolution.pmesh2domain import pmesh_to_domain_instance 
from pyvolution.pmesh2domain import pmesh_instance_to_domain_instance 
from caching import cache
import project
from create_mesh import create_mesh

#Data preparation
#Convert ASC 2 DEM 2 PTS using source data and store result in source data
demname = project.demname
demname = project.datadir + 'lakes_100' 
# just interested in the boundary/ mesh intereaction first off.
#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 CSIRO boundary
# NOTE: This function is dependent on a lot of files, and I haven't listed
# them, since there is so many and I don't think they will be changing.
cache(asc_csiro2sww,
      (project.bath_dir,
       project.elevation_dir,
       project.ucur_dir,
       project.vcur_dir,
       project.boundaryname + '.sww'),
       {'verbose': True,
       'minlat': -38.25,
       'maxlat': -37.7,
       'minlon': 147.0,
       'maxlon': 148.25,
       'mean_stage': project.tide,
       'zscale': 1,                 #Enhance tsunami
       'fail_on_NaN': False},
       #evaluate = True,
      verbose = True)


#Create the mesh ..

#cache this, with dependancy on create_mesh.py
mesh, triagle_count = create_mesh(10000000)


#meshname = project.meshname + '_%d.msh' %resolution
#meshname = project.meshdir + 'lakes_small_100.msh' 


#Setup domain

domain = cache(pmesh_instance_to_domain_instance, (mesh, Domain),
               #dependencies = [meshname],                     
               verbose = 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', project.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': Bd,'side': Bd, 'ocean': Bd} )  #Nothing


#Evolve
import time
t0 = time.time()

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)