source: trunk/anuga_validation/circular_island/simple_circular/circular.py @ 7877

Last change on this file since 7877 was 5460, checked in by kristy, 16 years ago
File size: 5.8 KB
Line 
1"""Simple water flow example using ANUGA
2
3Water flowing down a channel with more complex topography
4"""
5
6#------------------------------------------------------------------------------
7# Import necessary modules
8#------------------------------------------------------------------------------
9from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
10from anuga.shallow_water import Domain
11from anuga.shallow_water import Reflective_boundary
12from anuga.shallow_water import Dirichlet_boundary
13from anuga.shallow_water import Time_boundary
14from anuga.pmesh.mesh_interface import create_mesh_from_regions
15from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
16from math import tan, sqrt, sin, pi
17
18#------------------------------------------------------------------------------
19# Project file
20#------------------------------------------------------------------------------
21#Set up file structure
22anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
23meshes_dir = anuga_dir+'meshes'+sep
24meshname = 'circular_mesh.msh'
25output_dir = anuga_dir+'output'+sep
26
27
28##    #------------------------------------------------------------------------------
29##    # Copy scripts to time stamped output directory and capture screen
30##    # output to file
31##    #------------------------------------------------------------------------------
32##    print "Processor Name:",get_processor_name()
33##
34##    #copy script must be before screen_catcher
35##    #print kwargs
36##
37##    print 'output_dir', output_dir
38##    if myid == 0:
39##        copy_code_files(kwargs['output_dir'],__file__,
40##                 dirname(project.__file__)+sep+ project.__name__+'.py' )
41##
42##        store_parameters(**kwargs)
43##
44##    barrier()
45##
46##    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
47##
48##    print "Processor Name:",get_processor_name()
49##
50##    # filenames
51###    meshes_dir_name = project.meshes_dir_name+'.msh'
52##
53##    # creates copy of code in output dir
54##    print 'min triangles', project.trigs_min,
55##    print 'Note: This is generally about 20% less than the final amount'
56
57
58#------------------------------------------------------------------------------
59# Setup computational domain
60#------------------------------------------------------------------------------
61
62
63length = 30.
64width = 25.
65Cx = 12.96                  # centre of island on the x axis
66Cy = 13.8                   # centre of island on the y axis
67dx = dy = 1                 # Resolution: Length of subdivisions on both axes
68water_depth = 0.32          # Can be 0.32 or 0.42
69
70#boundary
71poly_domain = [[0,0],[length,0],[length,width],[0,width]]
72
73
74# exporting asc grid
75xmin = 0
76xmax = length   
77ymin = 0
78ymax = width
79
80#Create interior region
81Dome = [[(Cx)-4,(Cy)-4],[(Cx)+4,(Cy)-4,],
82        [(Cx)+4,(Cy)+4],[(Cx)-4,(Cy)+4]]
83
84remainder_res=1
85Dome_res = .01
86
87interior_dome = [[Dome, Dome_res]]
88
89#Create mesh
90
91print 'start create mesh from regions'
92 
93create_mesh_from_regions(poly_domain,
94                         boundary_tags={'wavemaker': [0], 'right': [1],
95                                        'top': [2], 'left': [3]},
96                         maximum_triangle_area = remainder_res,
97                         filename=meshname,
98                         interior_regions = interior_dome,
99                         use_cache=False,
100                         verbose=False)
101
102# Setup computational domain
103
104domain = Domain(meshes_dir+sep+meshname, use_cache=False, verbose = True)
105domain.set_name('circular')                  # Output name
106print 'memory usage before del domain',mem_usage()
107       
108    print domain.statistics()
109    print 'triangles',len(domain)
110
111
112#------------------------------------------------------------------------------
113# Setup initial conditions
114#------------------------------------------------------------------------------
115def topography(x,y):
116    """Complex topography defined by a function of vectors x and y
117    """
118   
119                                   
120    z= 0*x                      # defining z for all values other than the if statements
121    r= 3.6                      # radius, provided in document
122    angle = 14                  # angle, provided in document
123    h= r*tan(angle/57.2957795)  # finding height of cone if not truncated
124   
125
126    N = len(x)
127    for i in range(N):
128       
129        #truncated top
130        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <1.1**2:
131            z[i] += 0.625
132
133        # cone
134        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <r**2 and (x[i]-Cx)**2 + (y[i]-Cy)**2 >1.1**2:
135           z[i] = -(sqrt(((x[i]-Cx)**2+(y[i]-Cy)**2)/((r/h)**2))-h)
136    return z   
137
138domain.set_quantity('elevation', topography, verbose=True)  # Use function for elevation
139domain.set_quantity('friction', 0.01)         # Constant friction
140domain.set_quantity('stage',water_depth)   # Dry initial condition
141
142
143#------------------------------------------------------------------------------
144# Setup boundary conditions
145#------------------------------------------------------------------------------
146# Create boundary function from timeseries provided in file
147
148##boundary_filename="ts2cnew1_input_20_80sec_new"
149##prepare_timeboundary(boundary_filename+'.txt')
150##
151##function = file_function(boundary_filename+'.tms',
152##                         domain, verbose=True)
153def wave_form(t):
154    return 0.1*sin(2*pi*t/50.)
155
156# Create and assign boundary objects
157Bw = Dirichlet_boundary([water_depth, 0, 0])          #wall
158Bt = Transmissive_Momentum_Set_Stage_boundary(domain, wave_form) #wavemaker
159
160domain.set_boundary({'left': Bw, 'right': Bw, 'top': Bw, 'wavemaker': Bt})
161
162
163
164
165#------------------------------------------------------------------------------
166# Evolve system through time
167#------------------------------------------------------------------------------
168for t in domain.evolve(yieldstep = 0.2, finaltime = 100.0):
169    print domain.timestepping_statistics()
170
171
172   
Note: See TracBrowser for help on using the repository browser.