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

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 3.7 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.pyvolution.mesh_factory import rectangular
18from anuga.pyvolution.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.pyvolution.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 = 5,
36                          filename = 'island.msh' ,
37                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
38                          )
39
40
41
42#Create shallow water domain
43domain = Domain(mesh_filename = 'island.msh')
44domain.smooth = False
45domain.set_name('island')
46domain.set_default_order(2)
47
48
49#I tried to introduce this parameter top control the h-limiter,
50#but it doesn't remove the 'lapping water'
51#NEW (ON): I think it has been fixed (or at least reduced significantly) now
52#
53# beta_h == 1.0 means that the largest gradients (on h) are allowed
54# beta_h == 0.0 means that constant (1st order) gradients are introduced
55# on h. This is equivalent to the constant depth used previously.
56domain.beta_h = 0.9
57
58
59#------------------------------------------------------------------------------
60# Setup initial conditions
61#------------------------------------------------------------------------------
62
63def island(x, y):
64    z = 0*x
65    for i in range(len(x)):
66        z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
67
68        #z[i] += 0.5*exp( -((x[i]-10)**2 + (y[i]-10)**2)/50 )
69
70    return z
71
72def slump(x, y):
73    z = 0*x
74    for i in range(len(x)):
75        z[i] -= 0.7*exp( -((x[i]-10)**2 + (y[i]-10)**2)/200 )
76
77    return z
78
79#domain.set_quantity('friction', 0.1)  #Honky dory
80domain.set_quantity('friction', 100)     #Creep
81domain.set_quantity('elevation', island)
82domain.set_quantity('stage', 1)
83domain.max_timestep = 0.01
84
85
86#------------------------------------------------------------------------------
87# Setup boundary conditions (all reflective)
88#------------------------------------------------------------------------------
89
90Br = Reflective_boundary(domain)
91Bd = Dirichlet_boundary([1, 0, 0])
92
93domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
94domain.check_integrity()
95
96
97#------------------------------------------------------------------------------
98# Evolve system through time
99#------------------------------------------------------------------------------
100
101import time
102for t in domain.evolve(yieldstep = 1, finaltime = 100):
103    domain.write_time()
104    #if allclose(t, 100):
105    #    Q = domain.get_quantity('stage')
106    #    Q_s = Quantity(domain)
107    #    Q_s.set_values(slump)
108    #    domain.set_quantity('stage', Q + Q_s)
109    #print '    Volume: ', domain.get_quantity('stage').get_integral()
110
111print 'Done'
Note: See TracBrowser for help on using the repository browser.