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

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

Now test_parallel_sw_runup seems to be producing correct results (at
least on two processors)

File size: 10.1 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
35#--------------------------------------------------------------------------
36points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
37
38#if myid == 0:
39#    print 'points', points
40#    print 'vertices', vertices
41   
42domain = Domain(points, vertices, boundary) # Create domain
43
44#--------------------------------------------------------------------------
45# Setup initial conditions
46#--------------------------------------------------------------------------
47
48def topography(x,y): 
49    return -x/2                              # linear bed slope
50
51domain.set_quantity('elevation', topography) # Use function for elevation
52domain.set_quantity('friction', 0.0)         # Constant friction
53domain.set_quantity('stage', expression='elevation') # Dry initial stage
54
55
56#--------------------------------------------------------------------------
57# Create the parallel domain
58#--------------------------------------------------------------------------
59
60domain = distribute(domain, verbose=True)
61
62domain.set_name('runup')                    # Set sww filename
63domain.set_datadir('.')                     # Set output dir
64domain.set_maximum_allowed_speed(100)       #
65domain.set_quantities_to_be_stored(None)
66domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
67domain.H0 = 0 # Backwards compatibility (6/2/7)
68domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
69
70#------------------------------------------------------------------------------
71# Setup boundary conditions
72# This must currently happen *after* domain has been distributed
73#------------------------------------------------------------------------------
74
75Br = Reflective_boundary(domain)      # Solid reflective wall
76Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
77
78# Associate boundary tags with boundary objects
79domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
80
81
82
83#------------------------------------------------------------------------------
84# Evolve system through time
85#------------------------------------------------------------------------------
86
87interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
88gauge_values = []
89local_interpolation_points = []
90for i, point in enumerate(interpolation_points):
91    gauge_values.append([]) # Empty list for timeseries
92
93    if is_inside_polygon(point, domain.get_boundary_polygon()):
94
95        # FIXME: One point appears on multiple processes
96        # Need to get true boundary somehow
97       
98        #print 'P%d: point=[%f,%f]' %(myid, point[0], point[1])
99        local_interpolation_points.append(i)
100
101# Hack before we excluded ghosts.
102if numprocs == 2:
103    if myid == 0:
104        del local_interpolation_points[0]               
105        #local_interpolation_points = [1,2,3]
106if numprocs == 3:
107    if myid == 1:
108        del local_interpolation_points[0]
109if numprocs == 4:
110    if myid == 0:
111        del local_interpolation_points[1] #2
112        del local_interpolation_points[1] #3               
113    if myid == 3:
114        del local_interpolation_points[1]
115
116
117
118print 'P%d has points = %s' %(myid, local_interpolation_points)
119
120
121time = []
122
123for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
124    domain.write_time()
125   
126    # Record time series at known points
127    time.append(domain.get_time())
128   
129    stage = domain.get_quantity('stage')
130    w = stage.get_values(interpolation_points=interpolation_points)
131   
132    for i, _ in enumerate(interpolation_points):
133        gauge_values[i].append(w[i])
134
135if myid == 0:
136    for i, (x,y) in enumerate(interpolation_points):
137
138        try:
139            from pylab import *
140        except:
141            pass
142        else:
143            ion()
144            hold(False)
145            plot(time, gauge_values[i], 'r.-')
146            #time, predicted_gauge_values[i], 'k-')
147
148            title('Gauge %d (%f,%f)' %(i,x,y))
149            xlabel('time(s)')
150            ylabel('stage (m)')   
151            #legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
152            #savefig('Gauge_%d.png' %i, dpi = 300)
153
154        raw_input('Next')
155       
156
157
158# Reference from sequential version (also available as a
159# unit test in test_shallow_water_domain)
160# Added Friday 13 October 2006 by Ole
161
162        G0 = ensure_numeric([-0.20000000000000001, -0.20000000000000001, -0.19920600846161715, -0.19153647344085376, -0.19127622768281194, -0.1770671909675095, -0.16739412133181927, -0.16196038919122191, -0.15621633053131384, -0.15130021599977705, -0.13930978857215484, -0.19349274358263582, -0.19975307598803765, -0.19999897143103357, -0.1999999995532111, -0.19999999999949952, -0.19999999999949952, -0.19999999999949952, -0.19997270012494556, -0.19925805948554556, -0.19934513778450533, -0.19966484196394893, -0.1997352860102084, -0.19968260481750394, -0.19980280797303882, -0.19998804881822749, -0.19999999778075916, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167])
163
164        G1 = ensure_numeric([-0.29999999999999993, -0.29999588068034899, -0.29250047332330331, -0.28335081844518584, -0.26142206997410805, -0.22656028856329835, -0.21224087216745585, -0.19934324109114465, -0.1889857939783175, -0.18146311603911383, -0.17401078727434263, -0.15419361061257214, -0.16225060576782063, -0.19010941396999181, -0.20901161407004412, -0.21670683975774699, -0.21771386270738891, -0.21481284465869752, -0.21063120869004387, -0.20669243364582401, -0.20320707386714859, -0.19984087691926442, -0.19725417448019505, -0.19633783049069981, -0.19650494599999785, -0.19708316838336942, -0.19779309449413818, -0.19853070294429562, -0.19917342167307153, -0.19964814677795845, -0.19991627610824922, -0.20013162970144974, -0.20029864969405509, -0.20036259676501131, -0.20030682824965193, -0.20016105135750167, -0.19997664501985918, -0.19980185871568762, -0.19966836175417696, -0.19958856744312226, -0.19955954696194517, -0.19956950051110917, -0.19960377086336181, -0.19964885299433241, -0.19969427478531132, -0.19973301547655564, -0.19976121574277764, -0.19977765285688653, -0.19978315117522441, -0.19977994634841739, -0.19977101394878494])
165       
166        G2 = ensure_numeric([-0.40000000000000002, -0.39077401254732241, -0.33350466136630474, -0.29771023004255281, -0.27605439066140897, -0.25986156218997497, -0.24502185018573647, -0.231792624329521, -0.21981564668803993, -0.20870707082936543, -0.19877739883776599, -0.18980922837977957, -0.17308011674005838, -0.16306400164013773, -0.17798470933304333, -0.1929554075869116, -0.20236705191987037, -0.20695767560655007, -0.20841025876092567, -0.20792102174869989, -0.20655350005579293, -0.20492002526259828, -0.20310627026780645, -0.20105983335287836, -0.19937394565794653, -0.19853917506699659, -0.19836389977624452, -0.19850305023602796, -0.19877764028836831, -0.19910928131034669, -0.19943705712418805, -0.19970344172958865, -0.19991076989870474, -0.20010020127747646, -0.20025937787100062, -0.20035087292905965, -0.20035829921463297, -0.20029606557316171, -0.20019606915365515, -0.20009096093399206, -0.20000371608204368, -0.19994495432920584, -0.19991535665176338, -0.19990981826533513, -0.19992106419898723, -0.19994189853516578, -0.19996624091229293, -0.19998946016985167, -0.20000842303470234, -0.20002144460718174, -0.20002815561337187])
167       
168        G3 = ensure_numeric([-0.45000000000000001, -0.37631169657400332, -0.33000044342859486, -0.30586045469008522, -0.28843572253009941, -0.27215308978603808, -0.25712951540331219, -0.2431608296216613, -0.23032023651386374, -0.2184546873456619, -0.20735123704254332, -0.19740397194806389, -0.1859829564064375, -0.16675980728362105, -0.16951575032846536, -0.1832860872609344, -0.19485758939241243, -0.20231368291811427, -0.20625610376074754, -0.20758116241495619, -0.20721445402086161, -0.20603406830353785, -0.20450262808396991, -0.2026769581185151, -0.2007401212066364, -0.19931160535777592, -0.19863606301128725, -0.19848511940572691, -0.19860091042948352, -0.19885490669377764, -0.19916542732701112, -0.19946678238611959, -0.19971209594104697, -0.19991912886512292, -0.2001058430788881, -0.20024959409472989, -0.20032160254609382, -0.20031583165752354, -0.20025051539293123, -0.2001556115816068, -0.20005952955420872, -0.1999814429561611, -0.19992977821558131, -0.19990457708664208, -0.19990104785490476, -0.19991257153954825, -0.19993258231880562, -0.19995548502882532, -0.19997700760919687, -0.19999429663503748, -0.20000588800248761])
169
170       
171
172
173# Only compare those that belong to this process id
174G = [G0, G1, G2, G3]
175
176for i in local_interpolation_points:
177    msg = 'P%d, point #%d: Computed time series and reference time series are different: %s'\
178          %(myid, i, gauge_values[i]-G[i])
179    assert allclose(gauge_values[i], G[i]), msg
180
181print 'P%d completed succesfully using points = %s' %(myid, local_interpolation_points)
182
Note: See TracBrowser for help on using the repository browser.