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

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

Concentrating code

File size: 9.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
17import numpy as num
18
19from anuga.pmesh.mesh_interface import create_mesh_from_regions
20
21from anuga.utilities.numerical_tools import ensure_numeric
22from anuga.utilities.polygon import is_inside_polygon
23
24from anuga.interface import Domain
25from anuga.interface import Reflective_boundary
26from anuga.interface import Dirichlet_boundary
27from anuga.interface import Time_boundary
28from anuga.interface import Transmissive_boundary
29
30from anuga.interface import rectangular_cross
31
32from anuga_parallel.parallel_api import distribute, myid, numprocs
33
34
35#--------------------------------------------------------------------------
36# Setup computational domain on processor 0
37#--------------------------------------------------------------------------
38if myid == 0 :
39    points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
40
41    #print 'points', points
42    #print 'vertices', vertices
43   
44    domain = Domain(points, vertices, boundary) # Create domain
45else:
46    #So full domain is only constructed for the 0th processor (will save on memory)
47    domain = None
48
49#--------------------------------------------------------------------------
50# Setup initial conditions
51#--------------------------------------------------------------------------
52
53def topography(x,y): 
54    return -x/2                              # linear bed slope
55
56
57if myid == 0:
58    domain.set_quantity('elevation', topography) # Use function for elevation
59    domain.set_quantity('friction', 0.0)         # Constant friction
60    domain.set_quantity('stage', expression='elevation') # Dry initial stage
61
62
63#--------------------------------------------------------------------------
64# Create the parallel domain
65#--------------------------------------------------------------------------
66domain = distribute(domain, verbose=False)
67
68#--------------------------------------------------------------------------
69# Setup domain parameters
70#--------------------------------------------------------------------------
71domain.set_name('runup')                    # Set sww filename
72domain.set_datadir('.')                     # Set output dir
73
74domain.set_default_order(1)       
75domain.set_quantities_to_be_stored(None)
76
77
78#------------------------------------------------------------------------------
79# Setup boundary conditions
80# This must currently happen *AFTER* domain has been distributed
81#------------------------------------------------------------------------------
82
83Br = Reflective_boundary(domain)      # Solid reflective wall
84Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
85
86# Associate boundary tags with boundary objects
87domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
88
89
90
91#------------------------------------------------------------------------------
92# Evolve system through time
93#------------------------------------------------------------------------------
94
95interpolation_points = [[0.4,0.51], [0.6,0.51], [0.8,0.51], [0.9,0.51]]
96gauge_values = []
97tri_ids = []
98for i, point in enumerate(interpolation_points):
99    gauge_values.append([]) # Empty list for timeseries
100
101    #if is_inside_polygon(point, domain.get_boundary_polygon()):
102    try:
103        k = domain.get_triangle_containing_point(point)
104        print 'KKK',myid, k, domain.tri_full_flag[k]
105        if domain.tri_full_flag[k] == 1:
106            tri_ids.append(k)
107        else:
108            tri_ids.append(-1)           
109    except:
110        tri_ids.append(-2)
111
112print myid, tri_ids
113print myid, domain.tri_full_flag[tri_ids[0]]
114
115
116#print myid, domain.tri_full_flag
117
118
119
120print 'P%d has points = %s' %(myid, tri_ids)
121
122
123time = []
124
125for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
126    if myid == 0: domain.write_time()
127
128    # Record time series at known points
129    time.append(domain.get_time())
130   
131    stage = domain.get_quantity('stage')
132
133    for i in range(4):
134        if tri_ids[i] > -1:
135            gauge_values[i].append(stage.centroid_values[tri_ids[i]])
136
137
138
139
140G0 = [-0.19166666666666665, -0.19166666666666665, -0.19166666666666665, -0.18995135073536706, -0.17967620310033938, -0.17592302731546394, -0.17196188698352638, -0.16050146192256606, -0.15408782540491983, -0.15120360396251856, -0.15272836248974528, -0.18515612571281359, -0.19089313596572838, -0.19091930090793735, -0.19092564666264678, -0.19093961544356891, -0.19096113497288297, -0.19100903501132768, -0.19101957821753363, -0.19100676219446935, -0.19098186082883836, -0.19095020375973476, -0.19090160340144657, -0.19060035225639324, -0.19067236768901255, -0.19072092369650068, -0.19084725051273566, -0.19087147781405342, -0.19089622376237744, -0.19090954329991847, -0.19091808778328567, -0.19092562302462121, -0.19093120010640158, -0.19093455662311021, -0.1909359528144898, -0.19093587102253487, -0.19093485355514328, -0.19093340185413665, -0.19093189648362133, -0.19093071766205927, -0.19093005702136678, -0.19092994761677365, -0.19093031876104916, -0.19093104326472229, -0.19093197773630957, -0.19093299111521683, -0.19093398080658583, -0.19093487865230929, -0.19093564970049512, -0.1909362865237911, -0.19093680142438879]
141
142
143
144
145
146G1 = [-0.29166666666666669, -0.29166666666666669, -0.28685941525964553, -0.27057862587984083, -0.24118168337207896, -0.22173239931821093, -0.20976527589038554, -0.19971857867921899, -0.19173249780540819, -0.18487858041944857, -0.1780033786977564, -0.17215410368358283, -0.17155182647714193, -0.18655941293682041, -0.19977612728273281, -0.2068993409588821, -0.2095974217439926, -0.20961276275180968, -0.20804141273268489, -0.20583271698671743, -0.20358656706459732, -0.20159086143490923, -0.19997988026328414, -0.19903422512577046, -0.19862770982003813, -0.1984841836620021, -0.19852978191537585, -0.19875351136877187, -0.19912658475790856, -0.19955814714527378, -0.20002959678495849, -0.20032623461817603, -0.20045056162352595, -0.20047228698692546, -0.20043159564517379, -0.20035498975577987, -0.2002645931437364, -0.20017748770076391, -0.2001040624974407, -0.20005138031775113, -0.20002019418902675, -0.20000780866902876, -0.200010236268512, -0.20002255766703306, -0.2000398188461642, -0.20005785100286463, -0.2000737026875791, -0.20008570712104165, -0.20009331547005016, -0.20009681408763169, -0.20009700971531758]
147
148
149
150
151
152       
153G2 = [-0.39166666666666666, -0.38165411084088413, -0.33031114311718357, -0.29490845791170966, -0.27401187763459972, -0.25720432221568812, -0.24180370475026186, -0.22803625608700034, -0.21607642473560867, -0.2059173963924967, -0.19732372829977357, -0.18899569210124012, -0.18075177888497074, -0.1758458811698074, -0.18266970744342398, -0.19232195386649067, -0.199153652922028, -0.20292588217669899, -0.2046835952536466, -0.20515202619935408, -0.20476365635086183, -0.20384587744356614, -0.20266293067781077, -0.20143368989447288, -0.20039011693774927, -0.19970226213297601, -0.1992938050267305, -0.19906505599862009, -0.19899924116319651, -0.19908893561605648, -0.1993045458861577, -0.19960344595715143, -0.19989216726729134, -0.20008423737677972, -0.20018439752920908, -0.20022167063422625, -0.20021700109190652, -0.20018483837701212, -0.2001375242307058, -0.20008571798103356, -0.20003819123906194, -0.20000089191108147, -0.1999759332393749, -0.19996300730369299, -0.19996026140162992, -0.19996490481461601, -0.19997388431092483, -0.19998445575566867, -0.19999452242082996, -0.20000274678941998, -0.20000850192468142]
154
155
156
157       
158G3 = [-0.44166666666666665, -0.37566630156374725, -0.32907403981602412, -0.30514246507845882, -0.28736474524309508, -0.27062836126322631, -0.25490557254596713, -0.24015946060857596, -0.22672499664646506, -0.21493194009348415, -0.20488909335665775, -0.19594643601723874, -0.18708799192006037, -0.17881263621051871, -0.17829280002483108, -0.18600506642077144, -0.19389819408710718, -0.19941728321954535, -0.20269142693734038, -0.20431177320907121, -0.20475261565846536, -0.20438149942198236, -0.20350472280402784, -0.20238009732809753, -0.20123125943892334, -0.20030079794157421, -0.19967828320080519, -0.19928934067283594, -0.19907912341626033, -0.199034739569321, -0.19914305100399451, -0.19936563222000359, -0.19965218325198875, -0.19990375153648701, -0.20006953873503872, -0.20015884586614877, -0.20019243432639855, -0.20018670241799616, -0.20015514524359709, -0.2001098495335151, -0.2000610924794913, -0.20001718731119694, -0.19998313833853548, -0.19996070482302816, -0.1999494922032494, -0.1999476393706516, -0.1999524591476477, -0.19996109366774711, -0.19997102151684798, -0.19998034301901382, -0.1999878666221423]
159
160
161
162# Only compare those that belong to this process id
163G = [G0, G1, G2, G3]
164
165
166
167success = True
168
169for i in range(4):
170    if tri_ids[i] > -1:
171        #print 'myid = %g, allclose(gauge_values[%g], G%g) = %g' % (myid, i,i, allclose(gauge_values[i], G[i]))
172        success = success and num.allclose(gauge_values[i], G[i])
173
174
175
176if success:
177    print 'Successful completion on processor ',myid
178else:
179    print 'Failure on processor ',myid
Note: See TracBrowser for help on using the repository browser.