source: anuga_work/production/gippsland_2005/run_gippsland2.py @ 4294

#Convention for strings representing files
# #_file has the extention
#           #name does not have the extension

import time

from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
     dem2pts, asc_csiro2sww
from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance 
from caching import cache
from create_mesh2 import create_mesh
from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
     File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary
from anuga.pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 

# Import this last!  It stuffs up the loading of c extensions otherwise.
import project

#Data preparation
#Convert ASC 2 DEM 2 PTS using source data and store result in source data
# 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, project.demname[:-4], {'verbose': True},
      dependencies = [project.demname]      
      ,verbose = False
      )


# 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),
       {'verbose': True,
       'minlat': -38.25,
       'maxlat': -37.7,
       'minlon': 147.0,
       'maxlon': 148.25,
       'mean_stage': project.tide,
       'zscale': 1,                 #Enhance storm surge
       'fail_on_NaN': False}
      #,evaluate = True                  
      ,verbose = False
      )

meshname = project.meshname
pointname = project.pointname
mesh_elevname = project.mesh_elevname
outputname = project.outputname

t0 = time.time()
#Create the mesh without elevation data
# 100000 very course - 16,827 triangle mesh

meshname, triagle_count = cache(create_mesh,(1000000),
                                {'factor':1000000,
                                 'mesh_file':meshname,
                                 'triangles_in_name':True}
                                ,dependencies = ['create_mesh2.py']
                                #,evaluate = True      
                                )

mesh_elevname = mesh_elevname[:-9] + '_' + str(triagle_count) +  \
                mesh_elevname[-9:]
outputname = outputname[:-4] + '_' + str(triagle_count) + outputname[-4:]


#Add elevation data to the mesh
#fit_to_mesh_file(mesh_file, point_file, mesh_output_file, DEFAULT_ALPHA,
#                 verbose=True, expand_search=True, precrop=True)      
cache(fit_to_mesh_file,(meshname,
                 pointname,
                 mesh_elevname,
                 DEFAULT_ALPHA),
      {'verbose': True,
       'expand_search': True,
       'precrop': True}
      ,dependencies = [meshname, pointname]
      #,evaluate = True      
      ,verbose = False
      )
  
print 'Initialising the mesh took %.2f seconds' %(time.time()-t0)  

#Setup domain
domain = cache(pmesh_to_domain_instance, (mesh_elevname, Domain),
               dependencies = [mesh_elevname]                    
               ,verbose = False
               )               


domain.set_name(project.basename + '_%d' %triagle_count)
domain.set_datadir(project.outputdir)
domain.store = True
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.max_size = 3000
#Setup Boundary Conditions
print domain.get_boundary_tags()

print "****** run  ****"
print "project.boundaryname",project.boundaryname 
print "**********"

Bf = File_boundary(project.boundaryname, 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([project.tide,0,0])
#Bw = Time_boundary(domain=domain,
#                   f=lambda t: [(60<t<660)*4, 0, 0])

domain.set_boundary( {'back': Br,'side': Br, 'ocean': Bf} ) #CSIRO storm surge

#domain.set_boundary( {'back': Br,'side': Br, 'ocean': Br} ) #NO storm surge

#Evolve
t0 = time.time()

for t in domain.evolve(yieldstep = 900, finaltime = 432000):
#for t in domain.evolve(yieldstep = 10, finaltime = 30):
#                       skip_initial_step = True):
    domain.write_time()
    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')       

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