source: anuga_core/source/anuga_parallel/test_parallel_sw_runup.py @ 6721

Last change on this file since 6721 was 6721, checked in by steve, 15 years ago

Debugging parallel test function

File size: 9.4 KB
Line 
1#!/usr/bin/env python
2
3
4"""Simple water flow example using ANUGA
5
6Water driven up a linear slope and time varying boundary,
7similar to a beach environment
8
9This is a very simple test of the parallel algorithm using the simplified parallel API
10"""
11
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16
17from Numeric import allclose
18
19from anuga.pmesh.mesh_interface import create_mesh_from_regions
20from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
21from anuga.utilities.numerical_tools import ensure_numeric
22from anuga.utilities.polygon import is_inside_polygon
23
24from anuga.shallow_water import Domain
25from anuga.shallow_water import Reflective_boundary
26from anuga.shallow_water import Dirichlet_boundary
27from anuga.shallow_water import Time_boundary
28from anuga.shallow_water import Transmissive_boundary
29
30from parallel_api import distribute, myid, numprocs
31
32
33#--------------------------------------------------------------------------
34# Setup computational domain on processor 0
35#--------------------------------------------------------------------------
36if myid == 0 :
37    points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
38
39    #print 'points', points
40    #print 'vertices', vertices
41   
42    domain = Domain(points, vertices, boundary) # Create domain
43else:
44    #So full domain is only constructed for the 0th processor (will save on memory)
45    domain = None
46
47#--------------------------------------------------------------------------
48# Setup initial conditions
49#--------------------------------------------------------------------------
50
51def topography(x,y): 
52    return -x/2                              # linear bed slope
53
54
55if myid == 0:
56    domain.set_quantity('elevation', topography) # Use function for elevation
57    domain.set_quantity('friction', 0.0)         # Constant friction
58    domain.set_quantity('stage', expression='elevation') # Dry initial stage
59
60
61#--------------------------------------------------------------------------
62# Create the parallel domain
63#--------------------------------------------------------------------------
64domain = distribute(domain, verbose=False)
65
66#--------------------------------------------------------------------------
67# Setup domain parameters
68#--------------------------------------------------------------------------
69domain.set_name('runup')                    # Set sww filename
70domain.set_datadir('.')                     # Set output dir
71
72domain.set_default_order(1)       
73domain.set_quantities_to_be_stored(None)
74domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
75
76# FIXME (Ole): Need tests where this is commented out
77domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
78domain.H0 = 0 # Backwards compatibility (6/2/7)
79domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
80
81
82#------------------------------------------------------------------------------
83# Setup boundary conditions
84# This must currently happen *AFTER* domain has been distributed
85#------------------------------------------------------------------------------
86
87Br = Reflective_boundary(domain)      # Solid reflective wall
88Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
89
90# Associate boundary tags with boundary objects
91domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
92
93
94
95#------------------------------------------------------------------------------
96# Evolve system through time
97#------------------------------------------------------------------------------
98
99interpolation_points = [[0.4,0.51], [0.6,0.51], [0.8,0.51], [0.9,0.51]]
100gauge_values = []
101tri_ids = []
102for i, point in enumerate(interpolation_points):
103    gauge_values.append([]) # Empty list for timeseries
104
105    #if is_inside_polygon(point, domain.get_boundary_polygon()):
106    try:
107        k = domain.get_triangle_containing_point(point)
108        #print 'KKK',myid, k, domain.tri_full_flag[k]
109        if domain.tri_full_flag[k] == 1:
110            tri_ids.append(k)
111        else:
112            tri_ids.append(-1)           
113    except:
114        tri_ids.append(-2)
115
116#print myid, tri_ids
117#print myid, domain.tri_full_flag[tri_ids[0]]
118
119
120#print myid, domain.tri_full_flag
121
122
123
124print 'P%d has points = %s' %(myid, tri_ids)
125
126
127time = []
128
129for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
130    if myid == 0: domain.write_time()
131
132    # Record time series at known points
133    time.append(domain.get_time())
134   
135    stage = domain.get_quantity('stage')
136
137    for i in range(4):
138        if tri_ids[i] > -1:
139            gauge_values[i].append(stage.centroid_values[tri_ids[i]])
140
141
142
143
144G0 = [-0.19166666666666665, -0.19166666666666665, -0.1908726751282839, -0.18320314010752042, -0.18294289434947858, -0.17186692881708804, -0.16703039399924297, -0.16196038919122194, -0.15621633053131387, -0.15130021599977714, -0.13930978857215476, -0.18515941024930252, -0.19141974265470424, -0.19166563809769938, -0.19166666621987774, -0.19166666666616614, -0.19166666666616614, -0.19166666666616614, -0.19163936679820112, -0.19092472615522754, -0.19101180444286325, -0.19133150862916881, -0.1914019526842341, -0.19134927149082989, -0.19146947464147604, -0.19165471548476315, -0.19166666444742567, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831, -0.19166666666632831]
145
146G1 = [-0.29166666666666669, -0.29160604748111429, -0.28349212083663766, -0.26637122760395054, -0.25078182919408626, -0.22796566989512818, -0.21386588163016981, -0.20291676395026542, -0.19288397535979634, -0.18529863721918491, -0.17833464440180594, -0.16041289672597714, -0.15844675190030885, -0.18498595153509523, -0.20472977209518117, -0.21430915468432699, -0.21656349614594508, -0.2143710033505706, -0.21047592509949745, -0.20690277124897924, -0.20383621964071927, -0.2009312504622697, -0.19832906347760781, -0.19715218778669685, -0.19713511857310209, -0.19760258124272423, -0.19821523236790017, -0.19886735995941604, -0.19945089045017395, -0.19990124934465239, -0.20016164831872438, -0.20035891530554095, -0.20053280091253389, -0.20062879519403856, -0.20061453716771055, -0.20050516254700712, -0.20034457406881809, -0.2001804893235041, -0.20004787348381434, -0.19996313076460828, -0.19992645686226715, -0.19992808844809359, -0.19995462843450704, -0.19999336711919727, -0.20003430609120498, -0.20007059456815027, -0.20009823569348062, -0.20011560727571823, -0.20012297891084577, -0.20012200982886325, -0.2001151929510443]
147
148       
149G2 = [-0.39166666666666666, -0.38244067921398905, -0.33350466136630463, -0.29771023004255276, -0.27605439066140891, -0.25986156218997503, -0.24502185018573649, -0.23179262432952102, -0.21981564668803996, -0.20870707082936543, -0.19877739883776596, -0.18980922837977954, -0.1730801167400583, -0.16306400164063853, -0.17798470933316624, -0.19295540736943456, -0.20236705173335867, -0.20695767548229582, -0.20841025868691554, -0.20792102171307641, -0.20655350005113712, -0.2049200252815718, -0.20310627030300929, -0.20105983339290284, -0.19937394568595421, -0.1985391750807548, -0.19836389978207072, -0.19850305023816406, -0.19877764028800521, -0.19910928130800573, -0.19943705712074122, -0.19970344172627957, -0.19991076989513878, -0.20010020127344161, -0.20025937786728282, -0.20035087292654716, -0.20035829921368681, -0.20029606557358018, -0.2001960691549278, -0.20009096093557327, -0.20000371608352077, -0.19994495433034862, -0.1999153566524963, -0.19990981826568235, -0.19992106419902292, -0.19994189853498176, -0.19996624091198109, -0.19998946016949451, -0.20000842303436642, -0.20002144460691332, -0.20002815561319481]
150       
151G3 = [-0.44166666666666665, -0.37631169657400343, -0.33000044342859486, -0.30586045469008505, -0.28843572253009925, -0.27215308978603797, -0.25712951540331214, -0.2431608296216613, -0.23032023651386366, -0.21845468734566184, -0.20735123704254327, -0.19740397194806383, -0.18598295640643708, -0.16675980728412546, -0.1695157503295662, -0.18328608714846414, -0.19485758921965809, -0.2023136827680253, -0.20625610365424957, -0.20758116234856749, -0.20721445398906116, -0.20603406830101498, -0.20450262810357958, -0.2026769581530177, -0.20074012124164839, -0.19931160538111795, -0.1986360630235105, -0.19848511941161767, -0.1986009104317408, -0.19885490669335967, -0.19916542732465842, -0.19946678238299068, -0.19971209593781808, -0.19991912886144106, -0.20010584307496243, -0.20024959409137078, -0.20032160254398582, -0.2003158316568579, -0.20025051539345248, -0.20015561158283834, -0.20005952955567172, -0.19998144295750273, -0.19992977821660762, -0.19990457708729134, -0.19990104785520091, -0.19991257153955966, -0.19993258231861782, -0.19995548502852453, -0.19997700760886039, -0.19999429663472679, -0.20000588800224417]
152
153
154# Only compare those that belong to this process id
155G = [G0, G1, G2, G3]
156
157
158success = True
159
160for i in range(4):
161    if tri_ids[i] > -1:
162        #print 'myid = %g, allclose(gauge_values[%g], G%g) = %g' % (myid, i,i, allclose(gauge_values[i], G[i]))
163        success = success and allclose(gauge_values[i], G[i])
164
165
166if success:
167    print 'Successful completion on processor ',myid
168else:
169    print 'Failure on processor ',myid
Note: See TracBrowser for help on using the repository browser.