source: anuga_work/development/landslide_benchmark_problem/slide.py @ 5509

Last change on this file since 5509 was 5509, checked in by sexton, 16 years ago

landslide benchmark problem + data files

File size: 3.8 KB
Line 
1"""Simple water flow example using ANUGA
2
3Benchmark problem from teh third International Workshop on Long-Wave
4Runup Models
5http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=3
6
7"""
8
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
15from anuga.shallow_water import Domain
16from anuga.shallow_water import Reflective_boundary
17from anuga.shallow_water import Dirichlet_boundary
18from anuga.shallow_water import Time_boundary
19from anuga.shallow_water import Transmissive_boundary
20
21
22#------------------------------------------------------------------------------
23# Setup computational domain
24#------------------------------------------------------------------------------
25from anuga.pmesh.mesh_interface import create_mesh_from_regions
26bounding_polygon = [[0,10],[10,10],[10,0],[0,0]]
27remainder_res = 0.25
28meshname = 'test.msh'
29create_mesh_from_regions(bounding_polygon,
30                         boundary_tags={'top': [0],
31                                        'right': [1],
32                                        'bottom': [2],
33                                        'left': [3]},
34                         maximum_triangle_area=remainder_res,
35                         filename=meshname,
36                         #interior_regions=interior_regions,
37                         use_cache=False,
38                         verbose=True)
39
40from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
41from anuga.caching import cache
42
43domain = cache(pmesh_to_domain_instance,
44               (meshname, Domain),
45               dependencies = [meshname]) # Create domain
46domain.set_name('slide')                    # Output to file runup.sww
47domain.set_datadir('.')                     # Use current directory for output
48domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
49                                    'xmomentum',
50                                    'ymomentum'])   
51
52
53#------------------------------------------------------------------------------
54# Setup initial conditions
55#------------------------------------------------------------------------------
56
57from math import tan, radians
58beta = radians(5.7)
59delta = 1.
60mu = 0.01
61g = 9.81
62
63def topography(x,y):
64    return x*tan(beta)                       # H(x) = x tan(beta)
65
66domain.set_quantity('elevation', topography) # Use function for elevation
67domain.set_quantity('friction', 0.1)         # Constant friction
68domain.set_quantity('stage', 0.0)            # Constant negative initial stage
69
70from Numeric import sqrt, exp
71t = 0
72def ho(x,y):
73    return delta*exp(x)
74
75    #return delta*( -(2*sqrt(x*mu*mu/(delta*tan(beta))) - sqrt(g/delta)*mu*t)**2)
76#------------------------------------------------------------------------------
77# Setup boundary conditions
78#------------------------------------------------------------------------------
79
80from math import sin, pi, exp
81Br = Reflective_boundary(domain)      # Solid reflective wall
82Bt = Transmissive_boundary(domain)    # Continue all values on boundary
83Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
84Bw = Time_boundary(domain=domain,     # Time dependent boundary 
85                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])
86
87# Associate boundary tags with boundary objects
88domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
89
90
91#------------------------------------------------------------------------------
92# Evolve system through time
93#------------------------------------------------------------------------------
94
95for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0):
96    t1 = domain.get_time()
97    t = t1
98    domain.set_quantity('elevation', ho)
99    domain.write_time()
100   
101
102
Note: See TracBrowser for help on using the repository browser.