[2229] | 1 | """Validation of the AnuGA implementation of the shallow water wave equation. |
---|
| 2 | |
---|
| 3 | This script sets up LWRU2 benchmark with initial condition stated |
---|
| 4 | |
---|
| 5 | See also |
---|
| 6 | |
---|
| 7 | http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2 |
---|
| 8 | |
---|
| 9 | Depth at western boundary is d = 13.5 cm |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | |
---|
[2425] | 13 | use_variable_mesh = True #Use large variable mesh generated by create_mesh.py |
---|
| 14 | #use_variable_mesh = False #Use small structured mesh for speed |
---|
[2229] | 15 | |
---|
| 16 | |
---|
| 17 | import sys |
---|
| 18 | from os import sep |
---|
| 19 | sys.path.append('..'+sep+'..'+sep) #FIXME: Shouldn't be necessary |
---|
| 20 | |
---|
| 21 | def prepare_timeboundary(filename): |
---|
| 22 | """Converting benchmark 2 time series to NetCDF tms file. |
---|
| 23 | This is a 'throw-away' code taylor made for files like |
---|
| 24 | 'Benchmark_2_input.txt' from the LWRU2 benchmark |
---|
| 25 | """ |
---|
| 26 | |
---|
| 27 | print 'Preparing time boundary from %s' %filename |
---|
| 28 | from Numeric import array |
---|
| 29 | |
---|
| 30 | fid = open(filename) |
---|
| 31 | |
---|
| 32 | #Skip first line |
---|
| 33 | line = fid.readline() |
---|
| 34 | |
---|
| 35 | #Read remaining lines |
---|
| 36 | lines = fid.readlines() |
---|
| 37 | fid.close() |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | N = len(lines) |
---|
| 41 | T = zeros(N, Float) #Time |
---|
| 42 | Q = zeros(N, Float) #Values |
---|
| 43 | |
---|
| 44 | for i, line in enumerate(lines): |
---|
| 45 | fields = line.split() |
---|
| 46 | |
---|
| 47 | T[i] = float(fields[0]) |
---|
| 48 | Q[i] = float(fields[1]) |
---|
| 49 | |
---|
| 50 | |
---|
| 51 | #Create tms file |
---|
| 52 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 53 | |
---|
| 54 | outfile = filename[:-4] + '.tms' |
---|
| 55 | print 'Writing to', outfile |
---|
| 56 | fid = NetCDFFile(outfile, 'w') |
---|
| 57 | |
---|
| 58 | fid.institution = 'Geoscience Australia' |
---|
| 59 | fid.description = 'Input wave for Benchmark 2' |
---|
| 60 | fid.starttime = 0.0 |
---|
| 61 | fid.createDimension('number_of_timesteps', len(T)) |
---|
| 62 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
| 63 | fid.variables['time'][:] = T |
---|
| 64 | |
---|
| 65 | fid.createVariable('stage', Float, ('number_of_timesteps',)) |
---|
| 66 | fid.variables['stage'][:] = Q[:] |
---|
| 67 | |
---|
| 68 | fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) |
---|
| 69 | fid.variables['xmomentum'][:] = 0.0 |
---|
| 70 | |
---|
| 71 | fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) |
---|
| 72 | fid.variables['ymomentum'][:] = 0.0 |
---|
| 73 | |
---|
| 74 | fid.close() |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | # Module imports |
---|
| 78 | from pyvolution.shallow_water import Domain, Reflective_boundary,\ |
---|
| 79 | File_boundary, Transmissive_Momentum_Set_Stage_boundary |
---|
| 80 | from pyvolution.mesh_factory import rectangular_cross |
---|
| 81 | from pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
| 82 | from Numeric import array, zeros, Float, allclose |
---|
| 83 | import project |
---|
| 84 | from caching import cache |
---|
| 85 | |
---|
| 86 | #Preparing time boundary |
---|
| 87 | prepare_timeboundary(project.boundary_filename) |
---|
| 88 | |
---|
[2864] | 89 | #Preparing points (Only when using original Benchmark_2_Bathymetry.txt file) |
---|
| 90 | #from pyvolution.data_manager import xya2pts |
---|
| 91 | #xya2pts(project.bathymetry_filename, verbose = True, |
---|
| 92 | # z_func = lambda z: -z) |
---|
[2229] | 93 | |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | ####################### |
---|
| 97 | # Create Domain |
---|
| 98 | if use_variable_mesh is True: |
---|
| 99 | print 'Creating domain from', project.mesh_filename |
---|
| 100 | |
---|
[2866] | 101 | domain = Domain(project.mesh_filename, use_cache=True, verbose=True) |
---|
[2229] | 102 | |
---|
| 103 | else: |
---|
| 104 | print 'Creating regular mesh' |
---|
| 105 | N = 150 |
---|
| 106 | points, vertices, boundary = rectangular_cross(N, N/5*3, |
---|
| 107 | len1=5.448, len2=3.402) |
---|
| 108 | print 'Creating domain' |
---|
| 109 | |
---|
| 110 | #domain = Domain(points, vertices, boundary) |
---|
| 111 | domain = cache(Domain, (points, vertices, boundary)) |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | |
---|
[2535] | 115 | |
---|
[2229] | 116 | import sys, os |
---|
| 117 | base = os.path.basename(sys.argv[0]) |
---|
| 118 | domain.filename, _ = os.path.splitext(base) |
---|
| 119 | domain.default_order = 2 |
---|
| 120 | domain.store = True #Store for visualisation purposes |
---|
| 121 | |
---|
[2866] | 122 | #domain.check_integrity() |
---|
[2229] | 123 | print 'Number of triangles = ', len(domain) |
---|
| 124 | print 'The extent is ', domain.get_extent() |
---|
[2535] | 125 | print domain.statistics() |
---|
[2229] | 126 | |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | ####################### |
---|
| 130 | # Initial Conditions |
---|
| 131 | print 'Initial values' |
---|
| 132 | |
---|
| 133 | domain.set_quantity('elevation', |
---|
| 134 | filename = project.bathymetry_filename[:-4] + '.pts', |
---|
[2426] | 135 | alpha = 0.001, |
---|
[2229] | 136 | verbose = True, |
---|
| 137 | use_cache = True) |
---|
| 138 | |
---|
| 139 | domain.set_quantity('friction', 0.0) |
---|
| 140 | domain.set_quantity('stage', 0.0) |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | ###################### |
---|
| 145 | # Boundary Conditions |
---|
| 146 | # |
---|
| 147 | print 'Boundaries' |
---|
| 148 | Br = Reflective_boundary(domain) |
---|
| 149 | |
---|
| 150 | from pyvolution.util import file_function |
---|
| 151 | function = file_function(project.boundary_filename[:-4] + '.tms', |
---|
| 152 | domain, |
---|
| 153 | verbose = True) |
---|
| 154 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
| 155 | |
---|
| 156 | #Set boundary conditions |
---|
| 157 | if use_variable_mesh is True: |
---|
| 158 | domain.set_boundary({'wave': Bts, 'wall': Br}) |
---|
| 159 | else: |
---|
| 160 | domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br}) |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | |
---|
| 164 | ####################### |
---|
| 165 | # Evolve |
---|
| 166 | import time |
---|
| 167 | t0 = time.time() |
---|
| 168 | |
---|
| 169 | for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5): |
---|
| 170 | domain.write_time() |
---|
| 171 | |
---|
| 172 | print 'That took %.2f seconds' %(time.time()-t0) |
---|