[2229] | 1 | """ |
---|
| 2 | Author: John Jakeman |
---|
| 3 | Created: 12/12/2005 |
---|
| 4 | Main Cairns script using new interface |
---|
| 5 | """ |
---|
| 6 | |
---|
| 7 | #------------------------------- |
---|
| 8 | # Module imports |
---|
| 9 | #------------------------------- |
---|
| 10 | import sys, os |
---|
[3944] | 11 | from anuga.shallow_water import Domain, Reflective_boundary,\ |
---|
[2312] | 12 | File_boundary, Transmissive_Momentum_Set_Stage_boundary,\ |
---|
| 13 | Transmissive_boundary, Dirichlet_boundary |
---|
[3944] | 14 | #from anuga.pyvolution.mesh_factory import rectangular_cross |
---|
| 15 | #from anuga.pmesh2domain import pmesh_to_domain_instance |
---|
| 16 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[2229] | 17 | from Numeric import array, zeros, Float, allclose |
---|
[3434] | 18 | import cairns_project as project |
---|
| 19 | #import project |
---|
[2229] | 20 | from caching import cache |
---|
[3514] | 21 | from anuga.utilities.polygon import read_polygon, Polygon_function |
---|
[2229] | 22 | |
---|
[2312] | 23 | from conversion import convert_latlon_to_xycoords |
---|
| 24 | |
---|
[2229] | 25 | #------------------------------- |
---|
| 26 | # Domain |
---|
| 27 | #------------------------------- |
---|
| 28 | print 'Creating domain' |
---|
[2312] | 29 | |
---|
[2270] | 30 | """ |
---|
| 31 | M = 36 #float(len1)/m |
---|
| 32 | N = 10 #float(len2)/n |
---|
[2229] | 33 | |
---|
| 34 | points, elements, boundary = rectangular_cross(M, N, len1=6393693.9504678631, |
---|
| 35 | len2=1781111.2941870166, |
---|
| 36 | origin = (0.0, 0.0)) |
---|
| 37 | |
---|
| 38 | #points, elements, boundary = rectangular_cross(N, M, len1=58.0, len2=1.0, origin=(145.0,-24.0)) |
---|
| 39 | |
---|
[2270] | 40 | domain = cache(Domain, (points, elements, boundary)) |
---|
[2229] | 41 | |
---|
[2270] | 42 | """ |
---|
[2312] | 43 | |
---|
[3944] | 44 | #cache(pmesh_to_domain_instance, |
---|
| 45 | # (project.mesh_filename, Domain), |
---|
| 46 | # dependencies = [project.mesh_filename]) |
---|
| 47 | domain = Domain(project.mesh_filename, use_cache = True, verbose = True) |
---|
[2270] | 48 | |
---|
[2229] | 49 | domain.check_integrity() |
---|
| 50 | print 'Number of triangles = ', len(domain) |
---|
| 51 | print 'The extent is ', domain.get_extent() |
---|
| 52 | |
---|
| 53 | #--------------------------------------------------------- |
---|
| 54 | #Decide which quantities are to be stored at each timestep |
---|
| 55 | domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] |
---|
| 56 | |
---|
| 57 | #------------------------------- |
---|
| 58 | # Initial Conditions |
---|
| 59 | #------------------------------- |
---|
| 60 | print 'Initial values' |
---|
| 61 | |
---|
[2312] | 62 | latlong_origin = 145, -24 |
---|
| 63 | r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin) |
---|
| 64 | r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin) |
---|
| 65 | r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin) |
---|
| 66 | r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin) |
---|
| 67 | |
---|
| 68 | #p0 = read_polygon('cairns_fault.xya') |
---|
[2229] | 69 | #p0 = read_polygon('cairns_fault_degrees.xya') |
---|
[2312] | 70 | p0 = [r0, r1, r2, r3] |
---|
| 71 | domain.set_quantity('stage',Polygon_function([(p0,1.0)])) |
---|
[2229] | 72 | |
---|
| 73 | #print domain.quantities['stage'].vertex_values |
---|
| 74 | |
---|
| 75 | domain.set_quantity('elevation', |
---|
[2270] | 76 | filename = project.bathymetry_filename[:-4] + '.xya', |
---|
[2229] | 77 | alpha = 10.0, |
---|
| 78 | verbose = True, |
---|
| 79 | use_cache = True) |
---|
| 80 | |
---|
| 81 | # Need to increase elevation by the same amount as the stage |
---|
| 82 | # At present unsure how to achieve this |
---|
| 83 | #domain.set_quantity('elevation',Polygon_function([domain.quantities['elevation'].vertex_values+1000.0])) |
---|
| 84 | |
---|
| 85 | domain.set_quantity('friction', 0.0) |
---|
| 86 | |
---|
| 87 | #------------------------------- |
---|
| 88 | # Boundary conditions |
---|
| 89 | #------------------------------- |
---|
| 90 | print 'Boundaries' |
---|
| 91 | |
---|
| 92 | #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function) |
---|
[2312] | 93 | #Bd = Dirichlet_boundary([0,0,0]) |
---|
[2270] | 94 | Br = Reflective_boundary(domain) |
---|
| 95 | domain.set_boundary({'border': Br, 'border': Br, 'border': Br, 'border': Br}) |
---|
[2229] | 96 | |
---|
[2270] | 97 | #Bt = Transmissive_boundary(domain) |
---|
[2312] | 98 | #domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
[2229] | 99 | |
---|
| 100 | #------------------------------- |
---|
| 101 | # Setup domain runtime parameters |
---|
| 102 | #------------------------------- |
---|
[2270] | 103 | |
---|
[2312] | 104 | #domain.visualise = True |
---|
| 105 | #domain.visualise_color_stage = True |
---|
[2229] | 106 | |
---|
[2312] | 107 | #Set domain.visualise = False for large problem sizes |
---|
| 108 | domain.visualise = False |
---|
| 109 | domain.visualise_color_stage = False |
---|
[2270] | 110 | |
---|
[2229] | 111 | base = os.path.basename(sys.argv[0]) |
---|
[3846] | 112 | basename, _ = os.path.splitext(base) |
---|
| 113 | domain.set_name(basename) |
---|
[2229] | 114 | domain.default_order = 2 |
---|
| 115 | domain.store = True #Store for visualisation purposes |
---|
| 116 | |
---|
| 117 | #------------------------------- |
---|
| 118 | # Evolve |
---|
| 119 | #------------------------------- |
---|
| 120 | import time |
---|
| 121 | t0 = time.time() |
---|
[2270] | 122 | yieldstep = 60.0 |
---|
[2312] | 123 | finaltime = 21600.0 |
---|
[2229] | 124 | |
---|
| 125 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
| 126 | domain.write_time() |
---|
| 127 | |
---|
| 128 | print 'That took %.2f seconds' %(time.time()-t0) |
---|