""" Author: John Jakeman Created: 12/12/2005 Main Cairns script using new interface """ #------------------------------- # Module imports #------------------------------- import sys, os from pyvolution.shallow_water import Domain, Reflective_boundary,\ File_boundary, Transmissive_Momentum_Set_Stage_boundary,\ Transmissive_boundary, Dirichlet_boundary from pyvolution.mesh_factory import rectangular_cross from pyvolution.pmesh2domain import pmesh_to_domain_instance from Numeric import array, zeros, Float, allclose #import cairns_project import project from caching import cache from utilities.polygon import read_polygon, Polygon_function from conversion import convert_latlon_to_xycoords #------------------------------- # Domain #------------------------------- print 'Creating domain' """ M = 36 #float(len1)/m N = 10 #float(len2)/n points, elements, boundary = rectangular_cross(M, N, len1=6393693.9504678631, len2=1781111.2941870166, origin = (0.0, 0.0)) #points, elements, boundary = rectangular_cross(N, M, len1=58.0, len2=1.0, origin=(145.0,-24.0)) domain = cache(Domain, (points, elements, boundary)) """ domain = cache(pmesh_to_domain_instance, (project.mesh_filename, Domain), dependencies = [project.mesh_filename]) domain.check_integrity() print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() #--------------------------------------------------------- #Decide which quantities are to be stored at each timestep domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] #------------------------------- # Initial Conditions #------------------------------- print 'Initial values' latlong_origin = 145, -24 r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin) r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin) r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin) r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin) #p0 = read_polygon('cairns_fault.xya') #p0 = read_polygon('cairns_fault_degrees.xya') p0 = [r0, r1, r2, r3] domain.set_quantity('stage',Polygon_function([(p0,1.0)])) #print domain.quantities['stage'].vertex_values domain.set_quantity('elevation', filename = project.bathymetry_filename[:-4] + '.xya', alpha = 10.0, verbose = True, use_cache = True) # Need to increase elevation by the same amount as the stage # At present unsure how to achieve this #domain.set_quantity('elevation',Polygon_function([domain.quantities['elevation'].vertex_values+1000.0])) domain.set_quantity('friction', 0.0) #------------------------------- # Boundary conditions #------------------------------- print 'Boundaries' #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function) #Bd = Dirichlet_boundary([0,0,0]) Br = Reflective_boundary(domain) domain.set_boundary({'border': Br, 'border': Br, 'border': Br, 'border': Br}) #Bt = Transmissive_boundary(domain) #domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #------------------------------- # Setup domain runtime parameters #------------------------------- #domain.visualise = True #domain.visualise_color_stage = True #Set domain.visualise = False for large problem sizes domain.visualise = False domain.visualise_color_stage = False base = os.path.basename(sys.argv[0]) domain.filename, _ = os.path.splitext(base) domain.default_order = 2 domain.store = True #Store for visualisation purposes #------------------------------- # Evolve #------------------------------- import time t0 = time.time() yieldstep = 60.0 finaltime = 21600.0 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0)