source: inundation/examples/island.py @ 2731

Last change on this file since 2731 was 2731, checked in by steve, 18 years ago
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 pyvolution.mesh_factory import rectangular
18from pyvolution.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
19from pmesh.mesh_interface import create_mesh_from_regions
20from utilities.polygon import Polygon_function, read_polygon
21
22from 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.