1 | #!/usr/bin/env python |
---|
2 | # Test a run of the sequential shallow water domain against |
---|
3 | # a run of the parallel shallow water domain. |
---|
4 | # WARNING: This assumes that the command to run jobs is mpirun. |
---|
5 | # Tested with MPICH. |
---|
6 | |
---|
7 | #mesh_filename = "test-100.tsh" |
---|
8 | mesh_filename= "merimbula_10785_1.tsh" |
---|
9 | yieldstep = 1 |
---|
10 | finaltime = 90 |
---|
11 | quantity = 'stage' |
---|
12 | nprocs = 8 |
---|
13 | |
---|
14 | import unittest |
---|
15 | import os |
---|
16 | import sys |
---|
17 | import print_stats |
---|
18 | import pypar |
---|
19 | |
---|
20 | from Numeric import array, zeros, Float, take, nonzero |
---|
21 | from anuga.pyvolution.shallow_water import Domain |
---|
22 | from anuga.pyvolution.shallow_water import Reflective_boundary as sw_reflective_boundary |
---|
23 | from anuga.pyvolution.shallow_water import Transmissive_boundary as sw_transmissive_boundary |
---|
24 | from parallel_shallow_water import Parallel_Domain |
---|
25 | from parallel_shallow_water import Reflective_boundary as par_reflective_boundary |
---|
26 | from parallel_shallow_water import Transmissive_boundary as par_transmissive_boundary |
---|
27 | from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
28 | from anuga.utilities.norms import * |
---|
29 | from anuga.utilities.util_ext import double_precision |
---|
30 | from print_stats import print_test_stats, build_full_flag |
---|
31 | from pmesh_divide import pmesh_divide_metis |
---|
32 | from build_submesh import build_submesh |
---|
33 | from build_local import build_local_mesh |
---|
34 | from build_commun import send_submesh, rec_submesh, extract_hostmesh |
---|
35 | |
---|
36 | class Set_Stage: |
---|
37 | """Set an initial condition with constant water height, for x<x0 |
---|
38 | """ |
---|
39 | |
---|
40 | def __init__(self, x0=0.25, x1=0.5, h=1.0): |
---|
41 | self.x0 = x0 |
---|
42 | self.x1 = x1 |
---|
43 | self.h = h |
---|
44 | |
---|
45 | def __call__(self, x, y): |
---|
46 | return self.h*((x>self.x0)&(x<self.x1)) |
---|
47 | |
---|
48 | sys.path.append(".." + os.sep + "pyvolution") |
---|
49 | |
---|
50 | def parallel_test(): |
---|
51 | myid = pypar.rank() |
---|
52 | numprocs = pypar.size() |
---|
53 | proc_name = pypar.Get_processor_name() |
---|
54 | |
---|
55 | rect = zeros(4,Float) # Buffer for results |
---|
56 | if myid == 0: |
---|
57 | # Partition |
---|
58 | domain_full = pmesh_to_domain_instance(mesh_filename, Domain) |
---|
59 | rect = array(domain_full.xy_extent, Float) |
---|
60 | |
---|
61 | domain_full.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0)) |
---|
62 | #domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0)) |
---|
63 | domain_full.check_integrity() |
---|
64 | domain_full.smooth = False |
---|
65 | domain_full.reduction = min |
---|
66 | domain_full.checkpoint = False |
---|
67 | domain_full.visualise = False |
---|
68 | |
---|
69 | nodes, triangles, boundary, triangles_per_proc, quantities = \ |
---|
70 | pmesh_divide_metis(domain_full, numprocs) |
---|
71 | |
---|
72 | submesh = build_submesh(nodes, triangles, boundary,\ |
---|
73 | quantities, triangles_per_proc) |
---|
74 | |
---|
75 | for p in range(1, numprocs): |
---|
76 | send_submesh(submesh, triangles_per_proc, p) |
---|
77 | |
---|
78 | points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \ |
---|
79 | extract_hostmesh(submesh, triangles_per_proc) |
---|
80 | |
---|
81 | else: |
---|
82 | points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \ |
---|
83 | rec_submesh(0) |
---|
84 | |
---|
85 | pypar.broadcast(rect, 0) |
---|
86 | domain = Parallel_Domain(points, vertices, boundary, |
---|
87 | full_send_dict = full_send_dict, |
---|
88 | ghost_recv_dict = ghost_recv_dict) |
---|
89 | tri_full_flag = build_full_flag(domain, ghost_recv_dict) |
---|
90 | |
---|
91 | domain.default_order = 1 |
---|
92 | R = par_reflective_boundary(domain) |
---|
93 | domain.set_boundary({'outflow' : R, 'inflow' : R, 'inner' : R, 'exterior' : R, 'open' : R, 'ghost' : None}) |
---|
94 | domain.set_quantity('stage', quantities['stage']) |
---|
95 | domain.set_quantity('elevation', quantities['elevation']) |
---|
96 | domain.store = False |
---|
97 | |
---|
98 | l1list = [] |
---|
99 | l2list = [] |
---|
100 | linflist = [] |
---|
101 | l1norm = zeros(3, Float) |
---|
102 | l2norm = zeros(3, Float) |
---|
103 | linfnorm = zeros(3, Float) |
---|
104 | recv_norm = zeros(3, Float) |
---|
105 | # Evolution |
---|
106 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
107 | edges = take(domain.quantities[quantity].edge_values, nonzero(tri_full_flag)) |
---|
108 | l1norm[0] = l1_norm(edges[:, 0]) |
---|
109 | l1norm[1] = l1_norm(edges[:, 1]) |
---|
110 | l1norm[2] = l1_norm(edges[:, 2]) |
---|
111 | l2norm[0] = pow(l2_norm(edges[:,0]), 2) |
---|
112 | l2norm[1] = pow(l2_norm(edges[:,1]), 2) |
---|
113 | l2norm[2] = pow(l2_norm(edges[:,2]), 2) |
---|
114 | linfnorm[0] = linf_norm(edges[:,0]) |
---|
115 | linfnorm[1] = linf_norm(edges[:,1]) |
---|
116 | linfnorm[2] = linf_norm(edges[:,2]) |
---|
117 | if myid == 0: |
---|
118 | domain.write_time() |
---|
119 | for p in range(1, numprocs): |
---|
120 | pypar.receive(p, recv_norm) |
---|
121 | l1norm += recv_norm |
---|
122 | pypar.receive(p, recv_norm) |
---|
123 | l2norm += recv_norm |
---|
124 | pypar.receive(p, recv_norm) |
---|
125 | linfnorm[0] = max(linfnorm[0], recv_norm[0]) |
---|
126 | linfnorm[1] = max(linfnorm[1], recv_norm[1]) |
---|
127 | linfnorm[2] = max(linfnorm[2], recv_norm[2]) |
---|
128 | l1list.append(l1norm) |
---|
129 | l2norm[0] = pow(l2norm[0], 0.5) |
---|
130 | l2norm[1] = pow(l2norm[1], 0.5) |
---|
131 | l2norm[2] = pow(l2norm[2], 0.5) |
---|
132 | l2list.append(l2norm) |
---|
133 | linflist.append(linfnorm) |
---|
134 | else: |
---|
135 | pypar.send(l1norm, 0) |
---|
136 | pypar.send(l2norm, 0) |
---|
137 | pypar.send(linfnorm, 0) |
---|
138 | return (l1list, l2list, linflist) |
---|
139 | |
---|
140 | def sequential_test(): |
---|
141 | domain_full = pmesh_to_domain_instance(mesh_filename, Domain) |
---|
142 | |
---|
143 | domain_full.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0)) |
---|
144 | domain_full.check_integrity() |
---|
145 | domain_full.default_order = 1 |
---|
146 | domain_full.smooth = False |
---|
147 | domain_full.reduction = min |
---|
148 | domain_full.store = False |
---|
149 | domain_full.checkpoint = False |
---|
150 | domain_full.visualise = False |
---|
151 | R = sw_reflective_boundary(domain_full) |
---|
152 | domain_full.set_boundary({'outflow' : R, 'inflow' : R, 'inner' : R, 'exterior' : R, 'open' : R}) |
---|
153 | l1list = [] |
---|
154 | l2list = [] |
---|
155 | linflist = [] |
---|
156 | l1norm = zeros(3, Float) |
---|
157 | l2norm = zeros(3, Float) |
---|
158 | linfnorm = zeros(3, Float) |
---|
159 | for t in domain_full.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
160 | domain_full.write_time() |
---|
161 | edge = domain_full.quantities[quantity].edge_values |
---|
162 | l1norm[0] = l1_norm(edge[:,0]) |
---|
163 | l1norm[1] = l1_norm(edge[:,1]) |
---|
164 | l1norm[2] = l1_norm(edge[:,2]) |
---|
165 | l2norm[0] = l2_norm(edge[:,0]) |
---|
166 | l2norm[1] = l2_norm(edge[:,1]) |
---|
167 | l2norm[2] = l2_norm(edge[:,2]) |
---|
168 | linfnorm[0] = linf_norm(edge[:,0]) |
---|
169 | linfnorm[1] = linf_norm(edge[:,1]) |
---|
170 | linfnorm[2] = linf_norm(edge[:,2]) |
---|
171 | l1list.append(l1norm) |
---|
172 | l2list.append(l2norm) |
---|
173 | linflist.append(linfnorm) |
---|
174 | return (l1list, l2list, linflist) |
---|
175 | |
---|
176 | # Test an 8-way run of the shallow water equations |
---|
177 | # against the sequential code. |
---|
178 | |
---|
179 | class Test_Parallel_Sw(unittest.TestCase): |
---|
180 | def testParallelSw(self): |
---|
181 | print "Expect this test to fail if not run from the parallel directory." |
---|
182 | result = os.system("mpirun -np %d python test_parallel_sw.py" % nprocs) |
---|
183 | assert_(result == 0) |
---|
184 | |
---|
185 | # Because we are doing assertions outside of the TestCase class |
---|
186 | # the PyUnit defined assert_ function can't be used. |
---|
187 | def assert_(condition, msg="Assertion Failed"): |
---|
188 | if condition == False: |
---|
189 | pypar.finalize() |
---|
190 | raise AssertionError, msg |
---|
191 | |
---|
192 | if __name__=="__main__": |
---|
193 | if pypar.size() == 1: |
---|
194 | runner = unittest.TextTestRunner() |
---|
195 | suite = unittest.makeSuite(Test_Parallel_Sw, 'test') |
---|
196 | runner.run(suite) |
---|
197 | else: |
---|
198 | if pypar.rank() == 0: |
---|
199 | l1norm_seq, l2norm_seq, linfnorm_seq = sequential_test() |
---|
200 | l1norm_par, l2norm_par, linfnorm_par = parallel_test() |
---|
201 | if pypar.rank() == 0: |
---|
202 | assert_(len(l1norm_seq) == len(l1norm_par)) |
---|
203 | assert_(len(l2norm_seq) == len(l2norm_par)) |
---|
204 | assert_(len(linfnorm_seq) == len(linfnorm_par)) |
---|
205 | assert_(len(l1norm_seq) == len(l2norm_seq)) |
---|
206 | assert_(len(l2norm_seq) == len(linfnorm_seq)) |
---|
207 | # Anything smaller than tol we consider to be 0. |
---|
208 | # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15) |
---|
209 | # * 10 to cater for rounding error in computation. |
---|
210 | tol = pow(10, -1 * (double_precision() - 1)) |
---|
211 | for x in range(len(l1norm_seq)): |
---|
212 | for y in range(3): |
---|
213 | # Calculate relative difference in the norms |
---|
214 | assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol) |
---|
215 | assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol) |
---|
216 | assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol) |
---|
217 | if x > 0: |
---|
218 | # Verify that the quantity is being conserved across iterations. |
---|
219 | assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol) |
---|
220 | assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol) |
---|
221 | assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol) |
---|
222 | assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol) |
---|
223 | assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol) |
---|
224 | assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol) |
---|
225 | pypar.finalize() |
---|
226 | |
---|