[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 | |
---|
[3647] | 12 | # Module imports |
---|
| 13 | from Numeric import array, zeros, Float, allclose |
---|
[2229] | 14 | |
---|
[3647] | 15 | from anuga.shallow_water import Domain |
---|
| 16 | from anuga.shallow_water import Reflective_boundary |
---|
| 17 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
| 18 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
[2229] | 19 | |
---|
| 20 | import project |
---|
| 21 | |
---|
| 22 | |
---|
[3647] | 23 | #------------------------- |
---|
| 24 | # Create Domain |
---|
| 25 | #------------------------- |
---|
[2229] | 26 | |
---|
[3647] | 27 | use_variable_mesh = True #Use large variable mesh generated by create_mesh.py |
---|
| 28 | #use_variable_mesh = False #Use small structured mesh for speed |
---|
[2229] | 29 | |
---|
| 30 | if use_variable_mesh is True: |
---|
| 31 | print 'Creating domain from', project.mesh_filename |
---|
| 32 | |
---|
[2866] | 33 | domain = Domain(project.mesh_filename, use_cache=True, verbose=True) |
---|
[2229] | 34 | else: |
---|
[3705] | 35 | print 'Creating regular from regular mesh' |
---|
[2229] | 36 | N = 150 |
---|
| 37 | points, vertices, boundary = rectangular_cross(N, N/5*3, |
---|
| 38 | len1=5.448, len2=3.402) |
---|
[3705] | 39 | domain = Domain(points, vertices, boundary) |
---|
[2229] | 40 | |
---|
[3716] | 41 | domain.set_name('okushiri_old_limiters') |
---|
[3647] | 42 | domain.set_default_order(2) |
---|
[3657] | 43 | domain.set_minimum_storable_height(0.001) |
---|
[3716] | 44 | |
---|
| 45 | # Set old (pre Sep 2006) defaults for limiters |
---|
| 46 | domain.beta_w = 0.9 |
---|
| 47 | domain.beta_w_dry = 0.9 |
---|
| 48 | domain.beta_uh = 0.9 |
---|
| 49 | domain.beta_uh_dry = 0.9 |
---|
| 50 | domain.beta_vh = 0.9 |
---|
| 51 | domain.beta_vh_dry = 0.9 |
---|
| 52 | |
---|
[3705] | 53 | domain.check_integrity() |
---|
[2535] | 54 | print domain.statistics() |
---|
[2229] | 55 | |
---|
| 56 | |
---|
[3647] | 57 | #------------------------- |
---|
[2229] | 58 | # Initial Conditions |
---|
[3647] | 59 | #------------------------- |
---|
| 60 | domain.set_quantity('friction', 0.0) |
---|
| 61 | domain.set_quantity('stage', 0.0) |
---|
[2229] | 62 | domain.set_quantity('elevation', |
---|
[3705] | 63 | filename=project.bathymetry_filename[:-4] + '.pts', |
---|
| 64 | alpha=0.02, |
---|
| 65 | verbose=True, |
---|
| 66 | use_cache=True) |
---|
[2229] | 67 | |
---|
[3647] | 68 | #------------------------- |
---|
[2229] | 69 | # Boundary Conditions |
---|
[3647] | 70 | #------------------------- |
---|
[2229] | 71 | Br = Reflective_boundary(domain) |
---|
| 72 | |
---|
| 73 | function = file_function(project.boundary_filename[:-4] + '.tms', |
---|
| 74 | domain, |
---|
[3705] | 75 | verbose=True) |
---|
[3647] | 76 | |
---|
[2229] | 77 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
| 78 | |
---|
| 79 | if use_variable_mesh is True: |
---|
| 80 | domain.set_boundary({'wave': Bts, 'wall': Br}) |
---|
| 81 | else: |
---|
| 82 | domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br}) |
---|
| 83 | |
---|
| 84 | |
---|
[3647] | 85 | #------------------------- |
---|
| 86 | # Evolve through time |
---|
| 87 | #------------------------- |
---|
[2229] | 88 | import time |
---|
| 89 | t0 = time.time() |
---|
| 90 | |
---|
[3705] | 91 | w_max = 0 |
---|
[2229] | 92 | for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5): |
---|
| 93 | domain.write_time() |
---|
| 94 | |
---|
[3705] | 95 | w = domain.get_maximum_inundation_elevation() |
---|
| 96 | x, y = domain.get_maximum_inundation_location() |
---|
| 97 | print ' Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y) |
---|
| 98 | print |
---|
| 99 | |
---|
| 100 | if w > w_max: |
---|
| 101 | w_max = w |
---|
| 102 | x_max = x |
---|
| 103 | y_max = y |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | print '**********************************************' |
---|
| 107 | print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max) |
---|
| 108 | print '**********************************************' |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
[2229] | 112 | print 'That took %.2f seconds' %(time.time()-t0) |
---|