# source:anuga_work/development/demos/island.py@5440

Last change on this file since 5440 was 5440, checked in by ole, 16 years ago

Updated code investigating creep - there is still a little bit, but beta_h has no influence on this and can be retired as proposed in ticket:194

File size: 4.9 KB
Line
1"""Example of shallow water wave equation.
2
3Island surrounded by water.
4This example investigates onshore 'creep'
5
6Creep is defined as water moving uphill over many timesteps.
7The initial adjustment at time = 0 is not 'creep'
8
9"""
10
11
12#------------------------------------------------------------------------------
13# Import necessary modules
14#------------------------------------------------------------------------------
15
16# Standard modules
17from math import exp
18
19# Application specific imports
20from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
21from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
22from anuga.pmesh.mesh_interface import create_mesh_from_regions
24
25from anuga.abstract_2d_finite_volumes.quantity import Quantity
26from Numeric import allclose
27
28#------------------------------------------------------------------------------
29# Setup computational domain
30#------------------------------------------------------------------------------
31name = 'island'
32mesh_filename =  name + '.tsh'
33if True:
34    #Create basic mesh
35    create_mesh_from_regions( [[0,0], [100,0], [100,100], [0,100]],
36                              boundary_tags = {'bottom': [0],
37                                               'right': [1],
38                                               'top': [2],
39                                               'left': [3]},
40                              maximum_triangle_area = 2,
41                              filename =  mesh_filename,
42                              interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 10.0)]
43                              #interior_holes=[[[50,25], [70,25], [70,75], [50,75]]],
44                              )
45
46
47
48
49
50# Select precomputed meshes here. In changeset:4725 there is a huge difference
51# in timestepping characteristics depending on the mesh used.
52# The mesh generated on windows takes 10 times as many timesteps.
53
54#Create shallow water domain
55domain = Domain(mesh_filename = mesh_filename)
56domain.smooth = False #True
57domain.set_name(name)
58#domain.set_name('island_not_unique')
59domain.default_order = 2
60
61print domain.statistics()
62
63#I tried to introduce this parameter top control the h-limiter,
64#but it doesn't remove the 'lapping water'
65#NEW (ON): I think it has been fixed (or at least reduced significantly) now
66#
67# beta_h == 1.0 means that the largest gradients (on h) are allowed
68# beta_h == 0.0 means that constant (1st order) gradients are introduced
69# on h. This is equivalent to the constant depth used previously.
70#domain.beta_h     = 0.5
71#domain.beta_w_dry = 0.0
72#domain.alpha_balance = 10.0
73#domain.minimum_allowed_height = 1.0e-4
74#domain.maximum_allowed_speed = 0
75#domain.minimum_storable_height = 1.0e-4
76
77#------------------------
78# Combinations with creep
79#------------------------
80
81#domain.maximum_allowed_speed = 0
82#domain.beta_h = 0.0
83#domain.tight_slope_limiters = 0
84
85#domain.maximum_allowed_speed = 0
86#domain.beta_h = 0.0
87#domain.H0 = 0.01
88#domain.tight_slope_limiters = 1
89
90#---------------------------
91# Combinations with less creep
92#---------------------------
93
94#domain.maximum_allowed_speed = 1
95#domain.beta_h = 0.0
96#domain.tight_slope_limiters = 0
97
98domain.maximum_allowed_speed = 0
99domain.beta_h = 0.0
100domain.tight_slope_limiters = True
101
102
103
104#------------------------------------------------------------------------------
105# Setup initial conditions
106#------------------------------------------------------------------------------
107
108def island(x, y):
109    z = 0*x
110    for i in range(len(x)):
111        z[i] = 20*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
112
113        #z[i] += 0.5*exp( -((x[i]-10)**2 + (y[i]-10)**2)/50 )
114
115    return z
116
117def slump(x, y):
118    z = 0*x
119    for i in range(len(x)):
120        z[i] -= 0.7*exp( -((x[i]-10)**2 + (y[i]-10)**2)/200 )
121
122    return z
123
124stage_value = 10.0
125#domain.set_quantity('friction', 0.1)  #Honky dory
126domain.set_quantity('friction', 0.01)     #Creep
127domain.set_quantity('elevation', island)
128domain.set_quantity('stage', stage_value)
129domain.max_timestep = 0.01
130
131
132
133#------------------------------------------------------------------------------
134# Setup boundary conditions (all reflective)
135#------------------------------------------------------------------------------
136
137Br = Reflective_boundary(domain)
138Bd = Dirichlet_boundary([stage_value, 0, 0])
139
140domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd, 'exterior': Br})
141domain.check_integrity()
142
143#------------------------------------------------------------------------------
144# Evolve system through time
145#------------------------------------------------------------------------------
146
147import time
148for t in domain.evolve(yieldstep = 1, finaltime = 25):
149    domain.write_time()
150    #if allclose(t, 100):
151    #    Q = domain.get_quantity('stage')
152    #    Q_s = Quantity(domain)
153    #    Q_s.set_values(slump)
154    #    domain.set_quantity('stage', Q + Q_s)
155    #print '    Volume: ', domain.get_quantity('stage').get_integral()
156
157print 'Done'
Note: See TracBrowser for help on using the repository browser.