1 | """ |
---|
2 | Script for debugging potential problem with the boundary |
---|
3 | condition: transmissive_momentum_set_stage |
---|
4 | |
---|
5 | When the wave reflects back towards the boundary excessive values |
---|
6 | for momentum are generated leading to degenerate timesteps. |
---|
7 | |
---|
8 | Ole Nielsen and Duncan Gray, GA - 2008 |
---|
9 | """ |
---|
10 | |
---|
11 | |
---|
12 | #---------------------------------------------------------------------------- |
---|
13 | # Import necessary modules |
---|
14 | #---------------------------------------------------------------------------- |
---|
15 | |
---|
16 | # Related major packages |
---|
17 | from anuga.shallow_water import Domain |
---|
18 | from anuga.shallow_water import Reflective_boundary |
---|
19 | from anuga.shallow_water import Dirichlet_boundary |
---|
20 | from anuga.shallow_water import Time_boundary |
---|
21 | from anuga.shallow_water import File_boundary |
---|
22 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
23 | |
---|
24 | |
---|
25 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
26 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
27 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
28 | |
---|
29 | |
---|
30 | |
---|
31 | # Scenario specific imports |
---|
32 | import create_mesh |
---|
33 | |
---|
34 | |
---|
35 | def elevation_function(x,y): |
---|
36 | from Numeric import zeros, size, Float |
---|
37 | slope = 4 ## Bit of a magic Number |
---|
38 | |
---|
39 | z = zeros(size(x), Float) |
---|
40 | for i in range(len(x)): |
---|
41 | if x[i] < slope: |
---|
42 | z[i] = 0.0 |
---|
43 | else: |
---|
44 | z[i] = (x[i]-slope)*0.1 |
---|
45 | return z |
---|
46 | |
---|
47 | |
---|
48 | #------------------------------------------------------------------------- |
---|
49 | # Model parameters |
---|
50 | #------------------------------------------------------------------------- |
---|
51 | friction=0.01 |
---|
52 | |
---|
53 | basename = 'debug_boundary' |
---|
54 | |
---|
55 | boundary_file = 'boundary.tms' |
---|
56 | #boundary_file = 'boundary_mellow.tms' |
---|
57 | |
---|
58 | |
---|
59 | |
---|
60 | #------------------------------------------------------------------------- |
---|
61 | # Create the triangular mesh |
---|
62 | #------------------------------------------------------------------------- |
---|
63 | |
---|
64 | xright = 19.0 |
---|
65 | ybottom = 0 |
---|
66 | ytop = 0.45 |
---|
67 | xleft = 0.0 |
---|
68 | |
---|
69 | #Very simple mesh |
---|
70 | points, elements, boundary = rectangular_cross(200, 5, |
---|
71 | len1=xright, |
---|
72 | len2=ytop, |
---|
73 | origin = (0.0, 0.0)) |
---|
74 | |
---|
75 | |
---|
76 | |
---|
77 | #------------------------------------------------------------------------- |
---|
78 | # Setup computational domain |
---|
79 | #------------------------------------------------------------------------- |
---|
80 | #domain = Domain(mesh_filename, use_cache = False, verbose = True) |
---|
81 | domain = Domain(points, elements, boundary, verbose = True) |
---|
82 | print domain.statistics() |
---|
83 | |
---|
84 | domain.set_name(basename) |
---|
85 | domain.set_datadir('.') |
---|
86 | |
---|
87 | |
---|
88 | #------------------------------------------------------------------------- |
---|
89 | # Setup initial conditions |
---|
90 | #------------------------------------------------------------------------- |
---|
91 | |
---|
92 | domain.set_quantity('stage', 0.06) |
---|
93 | domain.set_quantity('friction', friction) |
---|
94 | domain.set_quantity('elevation', elevation_function) |
---|
95 | |
---|
96 | |
---|
97 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
98 | |
---|
99 | # Create boundary function from timeseries provided in file |
---|
100 | function = file_function(boundary_file, domain, verbose=True) |
---|
101 | |
---|
102 | |
---|
103 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
104 | Br = Reflective_boundary(domain) |
---|
105 | domain.set_boundary( {'top': Br, 'right': Br, 'bottom':Br, 'left': Bts} ) |
---|
106 | |
---|
107 | |
---|
108 | #------------------------------------------------------------------------- |
---|
109 | # Evolve system through time |
---|
110 | #------------------------------------------------------------------------- |
---|
111 | |
---|
112 | for t in domain.evolve(0.1, 15): |
---|
113 | print domain.timestepping_statistics(track_speeds=False) |
---|
114 | print domain.boundary_statistics(tags=['left']) |
---|
115 | |
---|
116 | print 'finished' |
---|
117 | |
---|