source: anuga_core/source/anuga/examples/island.py @ 4026

Last change on this file since 4026 was 4026, checked in by jack, 17 years ago

Moved the vpython visualiser to obsolete_code and cleared out things that call it from other code.

File size: 4.0 KB
Line 
1"""Example of shallow water wave equation.
2
3Island surrounded by water.
4This example investigates onshore 'creep'
5
6"""
7
8
9#------------------------------------------------------------------------------
10# Import necessary modules
11#------------------------------------------------------------------------------
12
13# Standard modules
14from math import exp
15
16# Application specific imports
17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
18from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
19from anuga.pmesh.mesh_interface import create_mesh_from_regions
20from anuga.utilities.polygon import Polygon_function, read_polygon
21
22from anuga.abstract_2d_finite_volumes.quantity import Quantity
23from Numeric import allclose
24
25#------------------------------------------------------------------------------
26# Setup computational domain
27#------------------------------------------------------------------------------
28
29#Create basic mesh
30create_mesh_from_regions( [[0,0], [100,0], [100,100], [0,100]],
31                          boundary_tags = {'bottom': [0],
32                                           'right': [1],
33                                           'top': [2],
34                                           'left': [3]},
35                          maximum_triangle_area = 1.0,
36                          filename = 'island.msh' ,
37                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 1.0)]
38                          #interior_holes=[[[50,25], [70,25], [70,75], [50,75]]],
39                          )
40
41
42
43#Create shallow water domain
44domain = Domain(mesh_filename = 'island.msh')
45domain.smooth = False
46domain.set_name('island')
47domain.default_order = 2
48
49
50#I tried to introduce this parameter top control the h-limiter,
51#but it doesn't remove the 'lapping water'
52#NEW (ON): I think it has been fixed (or at least reduced significantly) now
53#
54# beta_h == 1.0 means that the largest gradients (on h) are allowed
55# beta_h == 0.0 means that constant (1st order) gradients are introduced
56# on h. This is equivalent to the constant depth used previously.
57domain.beta_h     = 0.5
58domain.beta_w_dry = 0.0
59domain.alpha_balance = 10.0
60domain.minimum_allowed_height = 1.0e-4 
61domain.maximum_allowed_speed = 100.0 
62domain.minimum_storable_height = 1.0e-4 
63
64#------------------------------------------------------------------------------
65# Setup initial conditions
66#------------------------------------------------------------------------------
67
68def island(x, y):
69    z = 0*x
70    for i in range(len(x)):
71        z[i] = 20*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
72
73        #z[i] += 0.5*exp( -((x[i]-10)**2 + (y[i]-10)**2)/50 )
74
75    return z
76
77def slump(x, y):
78    z = 0*x
79    for i in range(len(x)):
80        z[i] -= 0.7*exp( -((x[i]-10)**2 + (y[i]-10)**2)/200 )
81
82    return z
83
84stage_value = 15.0
85#domain.set_quantity('friction', 0.1)  #Honky dory
86domain.set_quantity('friction', 0.01)     #Creep
87domain.set_quantity('elevation', island)
88domain.set_quantity('stage', stage_value)
89domain.max_timestep = 0.01
90
91
92#------------------------------------------------------------------------------
93# Setup boundary conditions (all reflective)
94#------------------------------------------------------------------------------
95
96Br = Reflective_boundary(domain)
97Bd = Dirichlet_boundary([stage_value, 0, 0])
98
99domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd, 'exterior': Br})
100domain.check_integrity()
101
102#------------------------------------------------------------------------------
103# Evolve system through time
104#------------------------------------------------------------------------------
105
106import time
107for t in domain.evolve(yieldstep = 1, finaltime = 100):
108    domain.write_time()
109    #if allclose(t, 100):
110    #    Q = domain.get_quantity('stage')
111    #    Q_s = Quantity(domain)
112    #    Q_s.set_values(slump)
113    #    domain.set_quantity('stage', Q + Q_s)
114    #print '    Volume: ', domain.get_quantity('stage').get_integral()
115
116print 'Done'
Note: See TracBrowser for help on using the repository browser.