1 | #!/usr/bin/env python |
---|
2 | |
---|
3 | """Test a run of the sequential shallow water domain against |
---|
4 | a run of the parallel shallow water domain. |
---|
5 | |
---|
6 | WARNING: This assumes that the command to run jobs is mpirun. |
---|
7 | Tested with MPICH and LAM (Ole) |
---|
8 | """ |
---|
9 | |
---|
10 | #------------------------------------------------------------------------------ |
---|
11 | # Import necessary modules |
---|
12 | #------------------------------------------------------------------------------ |
---|
13 | |
---|
14 | import unittest |
---|
15 | import os |
---|
16 | import sys |
---|
17 | import pypar |
---|
18 | |
---|
19 | from Numeric import allclose, array, zeros, Float, take, nonzero |
---|
20 | |
---|
21 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
22 | |
---|
23 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
24 | from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance |
---|
25 | |
---|
26 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
27 | from anuga.utilities.util_ext import double_precision |
---|
28 | from anuga.utilities.norms import l1_norm, l2_norm, linf_norm |
---|
29 | |
---|
30 | from anuga.shallow_water import Domain |
---|
31 | from anuga.shallow_water import Reflective_boundary |
---|
32 | from anuga.shallow_water import Dirichlet_boundary |
---|
33 | from anuga.shallow_water import Time_boundary |
---|
34 | from anuga.shallow_water import Transmissive_boundary |
---|
35 | |
---|
36 | |
---|
37 | from anuga_parallel.parallel_api import distribute, myid, numprocs |
---|
38 | |
---|
39 | |
---|
40 | #-------------------------------------------------------------------------- |
---|
41 | # Setup parameters |
---|
42 | #-------------------------------------------------------------------------- |
---|
43 | |
---|
44 | mesh_filename = "merimbula_10785_1.tsh" |
---|
45 | #mesh_filename = "test-100.tsh" |
---|
46 | yieldstep = 1 |
---|
47 | finaltime = 20 |
---|
48 | quantity = 'stage' |
---|
49 | nprocs = 4 |
---|
50 | |
---|
51 | #-------------------------------------------------------------------------- |
---|
52 | # Setup procedures |
---|
53 | #-------------------------------------------------------------------------- |
---|
54 | class Set_Stage: |
---|
55 | """Set an initial condition with constant water height, for x<x0 |
---|
56 | """ |
---|
57 | |
---|
58 | def __init__(self, x0=0.25, x1=0.5, h=1.0): |
---|
59 | self.x0 = x0 |
---|
60 | self.x1 = x1 |
---|
61 | self.h = h |
---|
62 | |
---|
63 | def __call__(self, x, y): |
---|
64 | return self.h*((x>self.x0)&(x<self.x1)) |
---|
65 | |
---|
66 | #-------------------------------------------------------------------------- |
---|
67 | # Setup test |
---|
68 | #-------------------------------------------------------------------------- |
---|
69 | def evolution_test(parallel=False): |
---|
70 | |
---|
71 | domain = pmesh_to_domain_instance(mesh_filename, Domain) |
---|
72 | domain.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0)) |
---|
73 | |
---|
74 | #-------------------------------------------------------------------------- |
---|
75 | # Create parallel domain if requested |
---|
76 | #-------------------------------------------------------------------------- |
---|
77 | |
---|
78 | if parallel: |
---|
79 | if myid == 0: print 'DISTRIBUTING PARALLEL DOMAIN' |
---|
80 | domain = distribute(domain, verbose=False) |
---|
81 | |
---|
82 | #------------------------------------------------------------------------------ |
---|
83 | # Setup boundary conditions |
---|
84 | # This must currently happen *after* domain has been distributed |
---|
85 | #------------------------------------------------------------------------------ |
---|
86 | domain.store = False |
---|
87 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
88 | |
---|
89 | domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br}) |
---|
90 | |
---|
91 | #------------------------------------------------------------------------------ |
---|
92 | # Setup diagnostic arrays |
---|
93 | #------------------------------------------------------------------------------ |
---|
94 | l1list = [] |
---|
95 | l2list = [] |
---|
96 | linflist = [] |
---|
97 | l1norm = zeros(3, Float) |
---|
98 | l2norm = zeros(3, Float) |
---|
99 | linfnorm = zeros(3, Float) |
---|
100 | recv_norm = zeros(3, Float) |
---|
101 | |
---|
102 | #------------------------------------------------------------------------------ |
---|
103 | # Evolution |
---|
104 | #------------------------------------------------------------------------------ |
---|
105 | if parallel: |
---|
106 | if myid == 0: print 'PARALLEL EVOLVE' |
---|
107 | else: |
---|
108 | print 'SEQUENTIAL EVOLVE' |
---|
109 | |
---|
110 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
111 | edges = take(domain.quantities[quantity].edge_values, nonzero(domain.tri_full_flag)) |
---|
112 | l1norm[0] = l1_norm(edges[:,0]) |
---|
113 | l1norm[1] = l1_norm(edges[:,1]) |
---|
114 | l1norm[2] = l1_norm(edges[:,2]) |
---|
115 | l2norm[0] = l2_norm(edges[:,0]) |
---|
116 | l2norm[1] = l2_norm(edges[:,1]) |
---|
117 | l2norm[2] = l2_norm(edges[:,2]) |
---|
118 | linfnorm[0] = linf_norm(edges[:,0]) |
---|
119 | linfnorm[1] = linf_norm(edges[:,1]) |
---|
120 | linfnorm[2] = linf_norm(edges[:,2]) |
---|
121 | if parallel: |
---|
122 | l2norm[0] = pow(l2norm[0], 2) |
---|
123 | l2norm[1] = pow(l2norm[1], 2) |
---|
124 | l2norm[2] = pow(l2norm[2], 2) |
---|
125 | if myid == 0: |
---|
126 | domain.write_time() |
---|
127 | |
---|
128 | #print edges[:,1] |
---|
129 | for p in range(1, numprocs): |
---|
130 | pypar.receive(p, recv_norm) |
---|
131 | l1norm += recv_norm |
---|
132 | pypar.receive(p, recv_norm) |
---|
133 | l2norm += recv_norm |
---|
134 | pypar.receive(p, recv_norm) |
---|
135 | linfnorm[0] = max(linfnorm[0], recv_norm[0]) |
---|
136 | linfnorm[1] = max(linfnorm[1], recv_norm[1]) |
---|
137 | linfnorm[2] = max(linfnorm[2], recv_norm[2]) |
---|
138 | |
---|
139 | l2norm[0] = pow(l2norm[0], 0.5) |
---|
140 | l2norm[1] = pow(l2norm[1], 0.5) |
---|
141 | l2norm[2] = pow(l2norm[2], 0.5) |
---|
142 | |
---|
143 | l1list.append(l1norm) |
---|
144 | l2list.append(l2norm) |
---|
145 | linflist.append(linfnorm) |
---|
146 | else: |
---|
147 | pypar.send(l1norm, 0) |
---|
148 | pypar.send(l2norm, 0) |
---|
149 | pypar.send(linfnorm, 0) |
---|
150 | else: |
---|
151 | domain.write_time() |
---|
152 | l1list.append(l1norm) |
---|
153 | l2list.append(l2norm) |
---|
154 | linflist.append(linfnorm) |
---|
155 | |
---|
156 | |
---|
157 | return (l1list, l2list, linflist) |
---|
158 | |
---|
159 | # Test an 8-way run of the shallow water equations |
---|
160 | # against the sequential code. |
---|
161 | |
---|
162 | class Test_Parallel_Sw(unittest.TestCase): |
---|
163 | def testParallelSw(self): |
---|
164 | print "Expect this test to fail if not run from the parallel directory." |
---|
165 | result = os.system("mpirun -np %d python test_parallel_sw.py" % nprocs) |
---|
166 | assert_(result == 0) |
---|
167 | |
---|
168 | # Because we are doing assertions outside of the TestCase class |
---|
169 | # the PyUnit defined assert_ function can't be used. |
---|
170 | def assert_(condition, msg="Assertion Failed"): |
---|
171 | if condition == False: |
---|
172 | #pypar.finalize() |
---|
173 | raise AssertionError, msg |
---|
174 | |
---|
175 | if __name__=="__main__": |
---|
176 | if numprocs == 1: |
---|
177 | runner = unittest.TextTestRunner() |
---|
178 | suite = unittest.makeSuite(Test_Parallel_Sw, 'test') |
---|
179 | runner.run(suite) |
---|
180 | else: |
---|
181 | |
---|
182 | pypar.barrier() |
---|
183 | if myid == 0: |
---|
184 | print 'SEQUENTIAL START' |
---|
185 | l1norm_seq, l2norm_seq, linfnorm_seq = evolution_test(parallel=False) |
---|
186 | |
---|
187 | pypar.barrier() |
---|
188 | if myid ==0: |
---|
189 | print 'PARALLEL START' |
---|
190 | |
---|
191 | l1norm_par, l2norm_par, linfnorm_par = evolution_test(parallel=True) |
---|
192 | |
---|
193 | if myid == 0: |
---|
194 | assert_(len(l1norm_seq) == len(l1norm_par)) |
---|
195 | assert_(len(l2norm_seq) == len(l2norm_par)) |
---|
196 | assert_(len(linfnorm_seq) == len(linfnorm_par)) |
---|
197 | assert_(len(l1norm_seq) == len(l2norm_seq)) |
---|
198 | assert_(len(l2norm_seq) == len(linfnorm_seq)) |
---|
199 | # Anything smaller than tol we consider to be 0. |
---|
200 | # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15) |
---|
201 | # * 10 to cater for rounding error in computation. |
---|
202 | tol = pow(10, -1 * (double_precision() - 1)) |
---|
203 | for x in range(len(l1norm_seq)): |
---|
204 | for y in range(3): |
---|
205 | # Calculate relative difference in the norms |
---|
206 | assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol) |
---|
207 | assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol) |
---|
208 | assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol) |
---|
209 | if x > 0: |
---|
210 | # Verify that the quantity is being conserved across iterations. |
---|
211 | assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol) |
---|
212 | assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol) |
---|
213 | assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol) |
---|
214 | assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol) |
---|
215 | assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol) |
---|
216 | assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol) |
---|
217 | |
---|
218 | print 'Parallel test OK' |
---|
219 | |
---|
220 | |
---|