source: anuga_work/development/analytical_solutions/Analytical_solution_contracting_channel.py @ 5029

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

Fixing up directory name

File size: 3.9 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a symmetrical converging frictionless channel.
3
4Specific methods pertaining to the 2D shallow water equation
5are imported from shallow_water
6for use with the generic finite volume framework
7
8   Copyright 2005
9   Christopher Zoppou, Stephen Roberts
10   ANU
11
12Specific methods pertaining to the 2D shallow water equation
13are imported from shallow_water
14for use with the generic finite volume framework
15
16Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
17numerical vector named conserved_quantities.
18"""
19
20#---------------
21# Module imports
22import sys
23from os import sep
24sys.path.append('..'+sep+'pyvolution')
25
26from shallow_water import Transmissive_boundary, Reflective_boundary, \
27     Dirichlet_boundary
28from shallow_water import Constant_height, Domain
29from pmesh2domain import pmesh_to_domain_instance
30from mesh_factory import contracting_channel_cross
31
32#-------
33# Domain from a file
34# filename = 'converging_channel_30846.tsh'
35# print 'Creating domain from', filename
36# domain = pmesh_to_domain_instance(filename, Domain)
37
38######################
39# Domain created within python
40#
41Total_length = 50
42W_upstream = 5.
43W_downstream = 2.5
44L_1 = 5.
45L_2 = 11
46L_3 = Total_length - L_1 - L_2
47n = 5
48m = 50
49
50points, elements, boundary = \
51     contracting_channel_cross(m, n, W_upstream, W_downstream, L_1, L_2, L_3)
52domain = Domain(points, elements, boundary)
53
54
55print 'Number of triangles = ', len(domain)
56
57#----------------
58# Order of scheme
59domain.default_order = 2
60domain.smooth = True
61
62#-------------------------------------
63# Provide file name for storing output
64domain.store = True     #Store for visualisation purposes
65domain.format = 'sww'   #Native netcdf visualisation format
66domain.filename = 'contracting_channel_second-order'
67
68
69#----------------------------------------------------------
70# Decide which quantities are to be stored at each timestep
71domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
72
73#------------------------------------------------
74# This is for Visual Python
75domain.visualise = True
76
77#------------------------------------------
78# Reduction operation for get_vertex_values
79#from anuga.pyvolution.util import mean
80#domain.reduction = mean
81
82#------------------------
83# Set boundary Conditions
84tags = {}
85tags['left']   = Dirichlet_boundary([0.2, 1.2, 0.0])
86tags['top']    = Reflective_boundary(domain)
87tags['bottom'] = Reflective_boundary(domain)
88tags['right']  = Transmissive_boundary(domain)
89domain.set_boundary(tags)
90
91#----------------------
92# Set initial condition
93domain.set_quantity('elevation', 0.0)
94domain.set_quantity('stage', 0.2)
95
96# Use the inscribed circle with safety factor of 0.9 to establish the time step
97# domain.set_to_inscribed_circle(safety_factor=0.9)
98
99#----------
100# Evolution
101import time
102t0 = time.time()
103for t in domain.evolve(yieldstep = 5.0, finaltime = 50.0):
104    domain.write_time()
105
106print 'That took %.2f seconds' %(time.time()-t0)
107
108N = domain.number_of_elements
109
110Stage = domain.quantities['stage']
111stage = Stage.centroid_values
112XY = domain.centroid_coordinates
113
114# Calculate average
115average_stage  = 0.0
116n_points = 0
117for n in range(N):
118    if XY[n,0] > 35.0:
119        average_stage = average_stage + stage[n]
120        n_points = n_points + 1
121
122average_stage = average_stage/n_points
123
124#Standard Deviation
125sigma = 0.0
126max_stage = -999999.
127min_stage = 999999
128for n in range(N):
129    if XY[n,0] > 35.0:
130        sigma = sigma + (average_stage - stage[n])**2
131        if stage[n] > max_stage:
132            max_stage = stage[n]
133        if stage[n] < min_stage:
134            min_stage = stage[n]
135
136import math
137sigma = math.sqrt(sigma/(n_points-1))
138
139print L_2, average_stage, sigma, max_stage, min_stage, n_points
140
141domain.initialise_visualiser()
142domain.visualiser.update_all()
Note: See TracBrowser for help on using the repository browser.