[488] | 1 | """Example of shallow water wave equation. |
---|
| 2 | |
---|
| 3 | This is called Netherlands because it shows a dam with a gap in it and |
---|
| 4 | stylised housed behind it and below the water surface. |
---|
| 5 | |
---|
| 6 | """ |
---|
| 7 | |
---|
| 8 | ###################### |
---|
| 9 | # Module imports |
---|
| 10 | # |
---|
| 11 | from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ |
---|
| 12 | Transmissive_boundary, Time_boundary, Constant_height |
---|
| 13 | |
---|
| 14 | from mesh_factory import from_polyfile, rectangular |
---|
| 15 | from Numeric import array |
---|
| 16 | from math import sqrt |
---|
[2941] | 17 | #from least_squares import Interpolation |
---|
| 18 | from fit_interpolate.interpolate import Interpolate |
---|
[488] | 19 | |
---|
| 20 | |
---|
| 21 | print 'Creating domain' |
---|
| 22 | #data_points, _, data_values = from_polyfile('cornell_room_medres') |
---|
| 23 | #points, triangles, values = from_polyfile('hires2') |
---|
| 24 | data_points, _, data_values = from_polyfile('hires2') |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | #Regrid onto numerically stable mesh |
---|
| 28 | # |
---|
| 29 | #Compute regular mesh based on resolution and extent of data |
---|
| 30 | data_points = array(data_points) |
---|
| 31 | pmax = max(data_points) |
---|
| 32 | pmin = min(data_points) |
---|
| 33 | |
---|
| 34 | M = len(data_points) |
---|
| 35 | |
---|
| 36 | N = int(0.8*sqrt(M)) |
---|
| 37 | |
---|
| 38 | #print N |
---|
| 39 | |
---|
| 40 | mesh_points, vertices, boundary = rectangular(N, N, |
---|
| 41 | len1=pmax[0]-pmin[0], |
---|
| 42 | len2=pmax[1]-pmin[1], |
---|
| 43 | origin = pmin) |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | #Compute smooth surface on new mesh based on values from old (regrid) |
---|
| 47 | print 'Interp' |
---|
[2941] | 48 | interp = Interpolate(mesh_points, vertices, alpha=10) |
---|
| 49 | mesh_values = interp.fit( data_points, data_values) # this has not been tested |
---|
[488] | 50 | print 'Len mesh values', len(mesh_values) |
---|
| 51 | print 'Len mesh points', len(mesh_points) |
---|
| 52 | |
---|
| 53 | |
---|
| 54 | #Create shallow water domain |
---|
| 55 | print 'Creating domain' |
---|
| 56 | domain = Domain(mesh_points, vertices) #, boundary) |
---|
| 57 | |
---|
| 58 | domain.check_integrity() |
---|
| 59 | domain.default_order = 2 |
---|
| 60 | domain.smooth = True |
---|
| 61 | domain.reduction = min #Looks a lot better on top of steep slopes |
---|
| 62 | |
---|
| 63 | print "Number of triangles = ", len(domain) |
---|
| 64 | |
---|
| 65 | domain.visualise = False |
---|
| 66 | domain.checkpoint = False |
---|
| 67 | domain.store = True #Store for visualisation purposes |
---|
| 68 | domain.format = 'sww' #Native netcdf visualisation format |
---|
| 69 | import sys, os |
---|
| 70 | root, ext = os.path.splitext(sys.argv[0]) |
---|
| 71 | if domain.smooth is True: |
---|
| 72 | s = 'smooth' |
---|
| 73 | else: |
---|
| 74 | s = 'nonsmooth' |
---|
| 75 | domain.filename = root + '_' + s |
---|
| 76 | |
---|
| 77 | #Set bed-slope and friction |
---|
| 78 | manning = 0.0 |
---|
| 79 | |
---|
| 80 | print 'Field values' |
---|
| 81 | domain.set_quantity('elevation', mesh_values) |
---|
| 82 | domain.set_quantity('friction', manning) |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | ###################### |
---|
| 86 | # Boundary conditions |
---|
| 87 | # |
---|
| 88 | print 'Boundaries' |
---|
| 89 | Br = Reflective_boundary(domain) |
---|
| 90 | domain.set_boundary({'exterior': Br}) |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | ###################### |
---|
| 95 | #Initial condition |
---|
| 96 | # |
---|
| 97 | print 'Initial condition' |
---|
| 98 | |
---|
| 99 | #Define water height as a lump in one corner |
---|
| 100 | def height(x, y): |
---|
| 101 | from Numeric import zeros, Float |
---|
| 102 | |
---|
| 103 | N = len(x) |
---|
| 104 | assert N == len(y) |
---|
| 105 | |
---|
| 106 | xmin = min(x); xmax = max(x) |
---|
| 107 | ymin = min(y); ymax = max(y) |
---|
| 108 | |
---|
| 109 | xrange = xmax - xmin |
---|
| 110 | yrange = ymax - ymin |
---|
| 111 | |
---|
| 112 | z = zeros(N, Float) |
---|
| 113 | for i in range(N): |
---|
| 114 | if x[i] <= xmin + 0.25*xrange and y[i] <= ymin + 0.25*yrange: |
---|
| 115 | z[i] = 300 |
---|
| 116 | |
---|
| 117 | return z |
---|
| 118 | |
---|
[773] | 119 | domain.set_quantity('stage', height) |
---|
[488] | 120 | |
---|
| 121 | E = domain.quantities['elevation'].vertex_values |
---|
[773] | 122 | L = domain.quantities['stage'].vertex_values |
---|
| 123 | domain.set_quantity('stage', E+L) |
---|
[488] | 124 | |
---|
| 125 | #Evolve |
---|
| 126 | for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0): |
---|
| 127 | domain.write_time() |
---|
| 128 | |
---|
| 129 | |
---|
| 130 | |
---|