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

Last change on this file since 3719 was 3719, checked in by steve, 17 years ago

Adding holes to create_mesh_from_regions

File size: 3.8 KB
RevLine 
[2638]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
[3560]17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
[3563]18from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
[3514]19from anuga.pmesh.mesh_interface import create_mesh_from_regions
20from anuga.utilities.polygon import Polygon_function, read_polygon
[2638]21
[3560]22from anuga.abstract_2d_finite_volumes.quantity import Quantity
[2638]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]},
[3719]35                          maximum_triangle_area = 1.0,
[2649]36                          filename = 'island.msh' ,
[3719]37                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 1.0)]
38                          #interior_holes=[[50,25], [70,25], [70,75], [50,75]],
[2638]39                          )
40
41
42
43#Create shallow water domain
44domain = Domain(mesh_filename = 'island.msh')
45domain.smooth = False
46domain.set_name('island')
47domain.set_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.
[3719]57domain.beta_h     = 0.2
58domain.beta_w_dry = 0.2
[2638]59
60
61#------------------------------------------------------------------------------
62# Setup initial conditions
63#------------------------------------------------------------------------------
64
65def island(x, y):
66    z = 0*x
67    for i in range(len(x)):
68        z[i] = 8*exp( -((x[i]-50)**2 + (y[i]-50)**2)/100 )
69
70        #z[i] += 0.5*exp( -((x[i]-10)**2 + (y[i]-10)**2)/50 )
71
72    return z
73
74def slump(x, y):
75    z = 0*x
76    for i in range(len(x)):
77        z[i] -= 0.7*exp( -((x[i]-10)**2 + (y[i]-10)**2)/200 )
78
79    return z
80
[2648]81#domain.set_quantity('friction', 0.1)  #Honky dory
[3719]82domain.set_quantity('friction', 100.0)     #Creep
[2638]83domain.set_quantity('elevation', island)
84domain.set_quantity('stage', 1)
[2648]85domain.max_timestep = 0.01
[2638]86
87
88#------------------------------------------------------------------------------
89# Setup boundary conditions (all reflective)
90#------------------------------------------------------------------------------
91
92Br = Reflective_boundary(domain)
93Bd = Dirichlet_boundary([1, 0, 0])
94
95domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
96domain.check_integrity()
97
[3719]98domain.initialise_visualiser()
[2638]99#------------------------------------------------------------------------------
100# Evolve system through time
101#------------------------------------------------------------------------------
102
103import time
104for t in domain.evolve(yieldstep = 1, finaltime = 100):
105    domain.write_time()
106    #if allclose(t, 100):
107    #    Q = domain.get_quantity('stage')
108    #    Q_s = Quantity(domain)
109    #    Q_s.set_values(slump)
110    #    domain.set_quantity('stage', Q + Q_s)
111    #print '    Volume: ', domain.get_quantity('stage').get_integral()
112
113print 'Done'
Note: See TracBrowser for help on using the repository browser.