source: anuga_validation/analytical solutions/Analytical_solution_contracting_channel.py @ 3846

Last change on this file since 3846 was 3846, checked in by ole, 17 years ago

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

File size: 3.8 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
22from anuga.shallow_water import Transmissive_boundary, Reflective_boundary, \
23     Dirichlet_boundary
24from anuga.shallow_water import Domain
25from pmesh2domain import pmesh_to_domain_instance
26from mesh_factory import contracting_channel_cross
27
28#-------
29# Domain from a file
30# filename = 'converging_channel_30846.tsh'
31# print 'Creating domain from', filename
32# domain = pmesh_to_domain_instance(filename, Domain)
33
34######################
35# Domain created within python
36#
37Total_length = 50
38W_upstream = 5.
39W_downstream = 2.5
40L_1 = 5.
41L_2 = 11
42L_3 = Total_length - L_1 - L_2
43n = 5
44m = 50
45
46points, elements, boundary = \
47     contracting_channel_cross(m, n, W_upstream, W_downstream, L_1, L_2, L_3)
48domain = Domain(points, elements, boundary)
49
50
51print 'Number of triangles = ', len(domain)
52
53#----------------
54# Order of scheme
55domain.default_order = 2
56domain.smooth = True
57
58#-------------------------------------
59# Provide file name for storing output
60domain.store = True     #Store for visualisation purposes
61domain.format = 'sww'   #Native netcdf visualisation format
62domain.set_name('contracting_channel_second-order')
63
64
65#----------------------------------------------------------
66# Decide which quantities are to be stored at each timestep
67domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
68
69#------------------------------------------------
70# This is for Visual Python
71domain.visualise = True
72
73#------------------------------------------
74# Reduction operation for get_vertex_values
75#from anuga.pyvolution.util import mean
76#domain.reduction = mean
77
78#------------------------
79# Set boundary Conditions
80tags = {}
81tags['left']   = Dirichlet_boundary([0.2, 1.2, 0.0])
82tags['top']    = Reflective_boundary(domain)
83tags['bottom'] = Reflective_boundary(domain)
84tags['right']  = Transmissive_boundary(domain)
85domain.set_boundary(tags)
86
87#----------------------
88# Set initial condition
89domain.set_quantity('elevation', 0.0)
90domain.set_quantity('stage', 0.2)
91
92# Use the inscribed circle with safety factor of 0.9 to establish the time step
93# domain.set_to_inscribed_circle(safety_factor=0.9)
94
95#----------
96# Evolution
97import time
98t0 = time.time()
99for t in domain.evolve(yieldstep = 5.0, finaltime = 50.0):
100    domain.write_time()
101
102print 'That took %.2f seconds' %(time.time()-t0)
103
104N = domain.number_of_elements
105
106Stage = domain.quantities['stage']
107stage = Stage.centroid_values
108XY = domain.centroid_coordinates
109
110# Calculate average
111average_stage  = 0.0
112n_points = 0
113for n in range(N):
114    if XY[n,0] > 35.0:
115        average_stage = average_stage + stage[n]
116        n_points = n_points + 1
117
118average_stage = average_stage/n_points
119
120#Standard Deviation
121sigma = 0.0
122max_stage = -999999.
123min_stage = 999999
124for n in range(N):
125    if XY[n,0] > 35.0:
126        sigma = sigma + (average_stage - stage[n])**2
127        if stage[n] > max_stage:
128            max_stage = stage[n]
129        if stage[n] < min_stage:
130            min_stage = stage[n]
131
132import math
133sigma = math.sqrt(sigma/(n_points-1))
134
135print L_2, average_stage, sigma, max_stage, min_stage, n_points
136
137domain.initialise_visualiser()
138domain.visualiser.update_all()
Note: See TracBrowser for help on using the repository browser.