1 | """Simple water flow example using ANUGA |
---|
2 | |
---|
3 | Water driven up a linear slope and time varying boundary, |
---|
4 | similar to a beach environment |
---|
5 | """ |
---|
6 | |
---|
7 | #------------------------------------------------------------------------------ |
---|
8 | # Import necessary modules |
---|
9 | #------------------------------------------------------------------------------ |
---|
10 | |
---|
11 | import sys |
---|
12 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
13 | from anuga.shallow_water import Domain |
---|
14 | from anuga.shallow_water import Reflective_boundary |
---|
15 | from anuga.shallow_water import Dirichlet_boundary |
---|
16 | from anuga.shallow_water import Time_boundary |
---|
17 | from anuga.shallow_water import Transmissive_boundary |
---|
18 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
19 | from anuga.geospatial_data.geospatial_data import * |
---|
20 | from math import cos |
---|
21 | |
---|
22 | #------------------------------------------------------------------------------ |
---|
23 | # Setup computational domain |
---|
24 | #------------------------------------------------------------------------------ |
---|
25 | dx = 1000. |
---|
26 | dy = dx |
---|
27 | L = 100000. |
---|
28 | W = dx |
---|
29 | |
---|
30 | # structured mesh |
---|
31 | points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2)) |
---|
32 | |
---|
33 | domain = Domain(points, vertices, boundary) |
---|
34 | |
---|
35 | domain.set_timestepping_method('euler') |
---|
36 | domain.set_default_order(2) |
---|
37 | domain.set_name('myexample9') |
---|
38 | domain.set_datadir('.') # Use current directory for output |
---|
39 | |
---|
40 | domain.beta_w = 1.0 |
---|
41 | domain.beta_w_dry = 0.2 |
---|
42 | domain.beta_uh = 1.0 |
---|
43 | domain.beta_uh_dry = 0.2 |
---|
44 | domain.beta_vh = 1.0 |
---|
45 | domain.beta_vh_dry = 0.2 |
---|
46 | domain.beta_h = 1.0 |
---|
47 | |
---|
48 | #------------------------------------------------------------------------------ |
---|
49 | # Setup initial conditions |
---|
50 | #------------------------------------------------------------------------------ |
---|
51 | #domain.set_quantity('elevation', topography) # Use function for elevation |
---|
52 | domain.set_quantity('elevation',-100) |
---|
53 | domain.set_quantity('friction', 0.00) |
---|
54 | domain.set_quantity('stage', 0.0) |
---|
55 | |
---|
56 | #----------------------------------------------------------------------------- |
---|
57 | # Setup boundary conditions |
---|
58 | #------------------------------------------------------------------------------ |
---|
59 | from math import sin, pi, exp |
---|
60 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
61 | Bt = Transmissive_boundary(domain) # Continue all values on boundary |
---|
62 | Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values |
---|
63 | amplitude = 1 |
---|
64 | #Bw = Transmissive_Momentum_Set_Stage_boundary(domain=domain, |
---|
65 | Bw = Time_boundary(domain=domain, # Time dependent boundary |
---|
66 | ## Sine wave |
---|
67 | f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0]) |
---|
68 | ## Sawtooth? |
---|
69 | # f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0]) |
---|
70 | ## Sharp rise, linear fall |
---|
71 | # f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0]) |
---|
72 | # f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0]) |
---|
73 | # f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0]) |
---|
74 | |
---|
75 | # Associate boundary tags with boundary objects |
---|
76 | domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br}) |
---|
77 | |
---|
78 | |
---|
79 | #=============================================================================== |
---|
80 | from anuga.visualiser import RealtimeVisualiser |
---|
81 | vis = RealtimeVisualiser(domain) |
---|
82 | #vis.render_quantity_height("elevation", zScale=1, offset = 5.0, dynamic=False) |
---|
83 | vis.render_quantity_height("stage", zScale =10000, dynamic=True) |
---|
84 | #vis.colour_height_quantity('stage', (lambda q:q['stage'], -1.0, 1.0)) |
---|
85 | vis.colour_height_quantity('stage', (1.0, 0.5, 0.5)) |
---|
86 | vis.start() |
---|
87 | #=============================================================================== |
---|
88 | |
---|
89 | |
---|
90 | #------------------------------------------------------------------------------ |
---|
91 | # Evolve system through time |
---|
92 | #------------------------------------------------------------------------------ |
---|
93 | |
---|
94 | for t in domain.evolve(yieldstep = 20.0, finaltime = 10*40*60.): |
---|
95 | domain.write_time() |
---|
96 | vis.update() |
---|
97 | |
---|
98 | vis.evolveFinished() |
---|
99 | |
---|