#Convention for strings representing files # #_file has the extention # #name does not have the extension import time from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\ dem2pts, asc_csiro2sww from pyvolution.pmesh2domain import pmesh_to_domain_instance from caching import cache from create_mesh2 import create_mesh from pyvolution.shallow_water import Domain, Reflective_boundary,\ File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary from 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,(5000000), {'factor':1000, '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) #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