# source:anuga_core/source/anuga_parallel/test_parallel_sw_runup.py@3928

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

Parallel domains now store only full triangles in sww files.
Still need to remove ghost nodes.

File size: 9.8 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
64domain.set_maximum_allowed_speed(100)       #
65
66
67#------------------------------------------------------------------------------
68# Setup boundary conditions
69# This must currently happen *after* domain has been distributed
70#------------------------------------------------------------------------------
71
72Br = Reflective_boundary(domain)      # Solid reflective wall
73Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
74
75# Associate boundary tags with boundary objects
76domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
77
78
79
80#------------------------------------------------------------------------------
81# Evolve system through time
82#------------------------------------------------------------------------------
83
84interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
85gauge_values = []
86local_interpolation_points = []
87for i, point in enumerate(interpolation_points):
88    gauge_values.append([]) # Empty list for timeseries
89
90    if is_inside_polygon(point, domain.get_boundary_polygon()):
91
92        # FIXME: One point appears on multiple processes
93        # Need to get true boundary somehow
94
95        #print 'P%d: point=[%f,%f]' %(myid, point[0], point[1])
96        local_interpolation_points.append(i)
97
98# Hack before we excluded ghosts.
99if numprocs == 2:
100    if myid == 0:
101        del local_interpolation_points[0]
102        #local_interpolation_points = [1,2,3]
103if numprocs == 3:
104    if myid == 1:
105        del local_interpolation_points[0]
106if numprocs == 4:
107    if myid == 0:
108        del local_interpolation_points[1] #2
109        del local_interpolation_points[1] #3
110    if myid == 3:
111        del local_interpolation_points[1]
112
113
114
115print 'P%d has points = %s' %(myid, local_interpolation_points)
116
117
118time = []
119
120for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
121    domain.write_time()
122
123    # Record time series at known points
124    time.append(domain.get_time())
125
126    stage = domain.get_quantity('stage')
127    w = stage.get_values(interpolation_points=interpolation_points)
128
129    for i, _ in enumerate(interpolation_points):
130        gauge_values[i].append(w[i])
131
132
133for i, (x,y) in enumerate(interpolation_points):
134
135    try:
136        from pylab import *
137    except:
138        pass
139    else:
140        ion()
141        hold(False)
142        plot(time, gauge_values[i], 'r.-')
143        #time, predicted_gauge_values[i], 'k-')
144
145        title('Gauge %d (%f,%f)' %(i,x,y))
146        xlabel('time(s)')
147        ylabel('stage (m)')
148        #legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
149        #savefig('Gauge_%d.png' %i, dpi = 300)
150
151        raw_input('Next')
152
153
154
155# Reference from sequential version (also available as a
156# unit test in test_shallow_water_domain)
157# Added Friday 13 October 2006 by Ole
158
159G0 = ensure_numeric([-0.20000000000000001, -0.19999681443389281, -0.1986192343695303, -0.19147413648863046, -0.19132688908678019, -0.17642317476621105, -0.167376262630034, -0.16192452887426961, -0.15609171725778803, -0.15127107084302249, -0.14048864340360018, -0.19296484125327093, -0.19997006390580363, -0.19999999999937063, -0.19999999999937063, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19999999999938772, -0.19974288463035494, -0.19951636867991712, -0.19966301435195755, -0.19981082259800226, -0.19978575003960128, -0.19992942471933109, -0.19999999931029933, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989, -0.19999999999906989])
160
161G1 = ensure_numeric([-0.29999999999999993, -0.29988962537199199, -0.29293904425532025, -0.28329367722887888, -0.25999146407696289, -0.22613875068011896, -0.21190705052094994, -0.19900707995208217, -0.18876305176191882, -0.18132447501091936, -0.17395459512711151, -0.15562414200985644, -0.16212999953643359, -0.18964422820514618, -0.20871181844346975, -0.21672207791083464, -0.21774940291862779, -0.21482868050219833, -0.21057786776704043, -0.20649663432591045, -0.20294932949211578, -0.19974459897911329, -0.19733648772704043, -0.19641404599824669, -0.19654095699184146, -0.19709942852191994, -0.19780873983410741, -0.19853259125123518, -0.19916495938961168, -0.19965391267799168, -0.19993539587158982, -0.2001383705551133, -0.20029344332295113, -0.20035349748150011, -0.20029886541561631, -0.20015541958920294, -0.19997273066429103, -0.19979879448668514, -0.19966016997024041, -0.19957558009501869, -0.19955725674938532, -0.19958083002853366, -0.19961752462568647, -0.19965296611330258, -0.19968998132634594, -0.19972532942208607, -0.19975372922008239, -0.19977196116929855, -0.19977951443660594, -0.19977792107284789, -0.19976991595502003])
162
163G2 = ensure_numeric([-0.40000000000000002, -0.39011996186687281, -0.33359026016903887, -0.29757449757405952, -0.27594124995715791, -0.25970211955309436, -0.24482929492054245, -0.23156757139219822, -0.21956485769139392, -0.20844522129026694, -0.19856327660654355, -0.18962303467030903, -0.17371085465024955, -0.16429840256208336, -0.17793711732368575, -0.19287799702389993, -0.20236271260796762, -0.20700727993623128, -0.20847704371373174, -0.20796895600687262, -0.20653398626186478, -0.20480656169870676, -0.20295863990994492, -0.20100199602968896, -0.19940642689498472, -0.19858371478015749, -0.19838672154605322, -0.19851093923669558, -0.19878191998909323, -0.19910827645394291, -0.19943514333832094, -0.19971231361970535, -0.19992429278849655, -0.20010744405928019, -0.20025927002359642, -0.20034751667523681, -0.20035504591467249, -0.20029401385620157, -0.20019492358237226, -0.20008934249434918, -0.19999808924091636, -0.19993869218976712, -0.19991589568150098, -0.19991815777945968, -0.19993012995477188, -0.19994576118144997, -0.19996497193815974, -0.19998586151236197, -0.20000487253824847, -0.20001903000364174, -0.20002698661385457])
164
165G3 = ensure_numeric([-0.45000000000000001, -0.37713945714588398, -0.33029565026933816, -0.30598209033945367, -0.28847101155177313, -0.27211191064563195, -0.25701544058818926, -0.24298945948410997, -0.23010402733784807, -0.21820351802867713, -0.20709938367218383, -0.19719881806182216, -0.18568281604361933, -0.16828653906676322, -0.16977310768235579, -0.1832707289594605, -0.19483524345250974, -0.20233480051649216, -0.20630757214159207, -0.20763927857964531, -0.20724458160595791, -0.20599191745446047, -0.20438329669495012, -0.20256105512496606, -0.20071269486729407, -0.19934403619901719, -0.19866860191898347, -0.19849975056296071, -0.19860870923007437, -0.19885838217851401, -0.19916422433758982, -0.19946861981642039, -0.19972267778871666, -0.19993013816258154, -0.20011063428833351, -0.20024891930311628, -0.20031882555219671, -0.20031326268593497, -0.20024881068472311, -0.20015443214902759, -0.20005669097631221, -0.19997542564643309, -0.19992564006223304, -0.19990746148869892, -0.19990923999172872, -0.19991956416813192, -0.19993484556273733, -0.1999538628054662, -0.19997381636620407, -0.19999130900268777, -0.20000388227457688])
166
167
168# Only compare those that belong to this process id
169G = [G0, G1, G2, G3]
170
171for i in local_interpolation_points:
172    msg = 'P%d, point #%d: Computed time series and reference time series are different: %s'\
173          %(myid, i, gauge_values[i]-G[i])
174    assert allclose(gauge_values[i], G[i]), msg
175
176print 'P%d completed succesfully using points = %s' %(myid, local_interpolation_points)
177
Note: See TracBrowser for help on using the repository browser.