source: anuga_work/development/demos/island.py @ 4720

Last change on this file since 4720 was 4720, checked in by ole, 17 years ago

Played with the 'creep' problem. I am still not sure what is going on.

File size: 4.6 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
23from anuga.utilities.polygon import Polygon_function, read_polygon
24
25from anuga.abstract_2d_finite_volumes.quantity import Quantity
26from Numeric import allclose
27
28#------------------------------------------------------------------------------
29# Setup computational domain
30#------------------------------------------------------------------------------
31
32#Create basic mesh
33create_mesh_from_regions( [[0,0], [100,0], [100,100], [0,100]],
34                          boundary_tags = {'bottom': [0],
35                                           'right': [1],
36                                           'top': [2],
37                                           'left': [3]},
38                          maximum_triangle_area = 2,
39                          filename = 'island.msh' ,
40                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 10.0)]
41                          #interior_holes=[[[50,25], [70,25], [70,75], [50,75]]],
42                          )
43
44
45
46#Create shallow water domain
47domain = Domain(mesh_filename = 'island.msh')
48domain.smooth = False #True
49domain.set_name('island_unique')
50#domain.set_name('island_not_unique')
51domain.default_order = 2
52
53
54#I tried to introduce this parameter top control the h-limiter,
55#but it doesn't remove the 'lapping water'
56#NEW (ON): I think it has been fixed (or at least reduced significantly) now
57#
58# beta_h == 1.0 means that the largest gradients (on h) are allowed
59# beta_h == 0.0 means that constant (1st order) gradients are introduced
60# on h. This is equivalent to the constant depth used previously.
61#domain.beta_h     = 0.5
62#domain.beta_w_dry = 0.0
63#domain.alpha_balance = 10.0
64#domain.minimum_allowed_height = 1.0e-4
65#domain.maximum_allowed_speed = 0
66#domain.minimum_storable_height = 1.0e-4
67
68#------------------------
69# Combinations with creep
70#------------------------
71
72#domain.maximum_allowed_speed = 0
73#domain.beta_h = 0.0
74#domain.tight_slope_limiters = 0
75
76#domain.maximum_allowed_speed = 0
77#domain.beta_h = 0.0
78#domain.H0 = 0.01
79#domain.tight_slope_limiters = 1
80
81#---------------------------
82# Combinations with less creep ??
83#---------------------------
84
85domain.maximum_allowed_speed = 1
86domain.beta_h = 0.0
87domain.tight_slope_limiters = 0
88
89#domain.maximum_allowed_speed = 10
90#domain.beta_h = 0.0
91#domain.H0 = 0.01
92#domain.tight_slope_limiters = 1
93
94
95
96
97
98#------------------------------------------------------------------------------
99# Setup initial conditions
100#------------------------------------------------------------------------------
101
102def island(x, y):
103    z = 0*x
104    for i in range(len(x)):
105        z[i] = 20*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
106
107        #z[i] += 0.5*exp( -((x[i]-10)**2 + (y[i]-10)**2)/50 )
108
109    return z
110
111def slump(x, y):
112    z = 0*x
113    for i in range(len(x)):
114        z[i] -= 0.7*exp( -((x[i]-10)**2 + (y[i]-10)**2)/200 )
115
116    return z
117
118stage_value = 10.0
119#domain.set_quantity('friction', 0.1)  #Honky dory
120domain.set_quantity('friction', 0.01)     #Creep
121domain.set_quantity('elevation', island)
122domain.set_quantity('stage', stage_value)
123domain.max_timestep = 0.01
124
125
126
127#------------------------------------------------------------------------------
128# Setup boundary conditions (all reflective)
129#------------------------------------------------------------------------------
130
131Br = Reflective_boundary(domain)
132Bd = Dirichlet_boundary([stage_value, 0, 0])
133
134domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd, 'exterior': Br})
135domain.check_integrity()
136
137#------------------------------------------------------------------------------
138# Evolve system through time
139#------------------------------------------------------------------------------
140
141import time
142for t in domain.evolve(yieldstep = 1, finaltime = 25):
143    domain.write_time()
144    #if allclose(t, 100):
145    #    Q = domain.get_quantity('stage')
146    #    Q_s = Quantity(domain)
147    #    Q_s.set_values(slump)
148    #    domain.set_quantity('stage', Q + Q_s)
149    #print '    Volume: ', domain.get_quantity('stage').get_integral()
150
151print 'Done'
Note: See TracBrowser for help on using the repository browser.