| 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 shallow_water import Domain |
|---|
| 22 | from shallow_water import Reflective_boundary as sw_reflective_boundary |
|---|
| 23 | from 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 pmesh2domain import pmesh_to_domain_instance |
|---|
| 28 | from utilities.norms import * |
|---|
| 29 | from 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, extract_hostmesh |
|---|
| 33 | from build_local import build_local_mesh |
|---|
| 34 | from build_commun import send_submesh, rec_submesh |
|---|
| 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 | hostmesh = extract_hostmesh(submesh) |
|---|
| 79 | points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \ |
|---|
| 80 | build_local_mesh(hostmesh, 0, triangles_per_proc[0], numprocs) |
|---|
| 81 | |
|---|
| 82 | else: |
|---|
| 83 | points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \ |
|---|
| 84 | rec_submesh(0) |
|---|
| 85 | |
|---|
| 86 | pypar.broadcast(rect, 0) |
|---|
| 87 | domain = Parallel_Domain(points, vertices, boundary, |
|---|
| 88 | full_send_dict = full_send_dict, |
|---|
| 89 | ghost_recv_dict = ghost_recv_dict) |
|---|
| 90 | tri_full_flag = build_full_flag(domain, ghost_recv_dict) |
|---|
| 91 | |
|---|
| 92 | domain.default_order = 1 |
|---|
| 93 | R = par_reflective_boundary(domain) |
|---|
| 94 | domain.set_boundary({'outflow' : R, 'inflow' : R, 'inner' : R, 'exterior' : R, 'open' : R, 'ghost' : None}) |
|---|
| 95 | domain.set_quantity('stage', quantities['stage']) |
|---|
| 96 | domain.set_quantity('elevation', quantities['elevation']) |
|---|
| 97 | domain.store = False |
|---|
| 98 | |
|---|
| 99 | l1list = [] |
|---|
| 100 | l2list = [] |
|---|
| 101 | linflist = [] |
|---|
| 102 | l1norm = zeros(3, Float) |
|---|
| 103 | l2norm = zeros(3, Float) |
|---|
| 104 | linfnorm = zeros(3, Float) |
|---|
| 105 | recv_norm = zeros(3, Float) |
|---|
| 106 | # Evolution |
|---|
| 107 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
|---|
| 108 | edges = take(domain.quantities[quantity].edge_values, nonzero(tri_full_flag)) |
|---|
| 109 | l1norm[0] = l1_norm(edges[:, 0]) |
|---|
| 110 | l1norm[1] = l1_norm(edges[:, 1]) |
|---|
| 111 | l1norm[2] = l1_norm(edges[:, 2]) |
|---|
| 112 | l2norm[0] = pow(l2_norm(edges[:,0]), 2) |
|---|
| 113 | l2norm[1] = pow(l2_norm(edges[:,1]), 2) |
|---|
| 114 | l2norm[2] = pow(l2_norm(edges[:,2]), 2) |
|---|
| 115 | linfnorm[0] = linf_norm(edges[:,0]) |
|---|
| 116 | linfnorm[1] = linf_norm(edges[:,1]) |
|---|
| 117 | linfnorm[2] = linf_norm(edges[:,2]) |
|---|
| 118 | if myid == 0: |
|---|
| 119 | domain.write_time() |
|---|
| 120 | for p in range(1, numprocs): |
|---|
| 121 | pypar.receive(p, recv_norm) |
|---|
| 122 | l1norm += recv_norm |
|---|
| 123 | pypar.receive(p, recv_norm) |
|---|
| 124 | l2norm += recv_norm |
|---|
| 125 | pypar.receive(p, recv_norm) |
|---|
| 126 | linfnorm[0] = max(linfnorm[0], recv_norm[0]) |
|---|
| 127 | linfnorm[1] = max(linfnorm[1], recv_norm[1]) |
|---|
| 128 | linfnorm[2] = max(linfnorm[2], recv_norm[2]) |
|---|
| 129 | l1list.append(l1norm) |
|---|
| 130 | l2norm[0] = pow(l2norm[0], 0.5) |
|---|
| 131 | l2norm[1] = pow(l2norm[1], 0.5) |
|---|
| 132 | l2norm[2] = pow(l2norm[2], 0.5) |
|---|
| 133 | l2list.append(l2norm) |
|---|
| 134 | linflist.append(linfnorm) |
|---|
| 135 | else: |
|---|
| 136 | pypar.send(l1norm, 0) |
|---|
| 137 | pypar.send(l2norm, 0) |
|---|
| 138 | pypar.send(linfnorm, 0) |
|---|
| 139 | return (l1list, l2list, linflist) |
|---|
| 140 | |
|---|
| 141 | def sequential_test(): |
|---|
| 142 | domain_full = pmesh_to_domain_instance(mesh_filename, Domain) |
|---|
| 143 | |
|---|
| 144 | domain_full.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0)) |
|---|
| 145 | domain_full.check_integrity() |
|---|
| 146 | domain_full.default_order = 1 |
|---|
| 147 | domain_full.smooth = False |
|---|
| 148 | domain_full.reduction = min |
|---|
| 149 | domain_full.store = False |
|---|
| 150 | domain_full.checkpoint = False |
|---|
| 151 | domain_full.visualise = False |
|---|
| 152 | R = sw_reflective_boundary(domain_full) |
|---|
| 153 | domain_full.set_boundary({'outflow' : R, 'inflow' : R, 'inner' : R, 'exterior' : R, 'open' : R}) |
|---|
| 154 | l1list = [] |
|---|
| 155 | l2list = [] |
|---|
| 156 | linflist = [] |
|---|
| 157 | l1norm = zeros(3, Float) |
|---|
| 158 | l2norm = zeros(3, Float) |
|---|
| 159 | linfnorm = zeros(3, Float) |
|---|
| 160 | for t in domain_full.evolve(yieldstep = yieldstep, finaltime = finaltime): |
|---|
| 161 | domain_full.write_time() |
|---|
| 162 | edge = domain_full.quantities[quantity].edge_values |
|---|
| 163 | l1norm[0] = l1_norm(edge[:,0]) |
|---|
| 164 | l1norm[1] = l1_norm(edge[:,1]) |
|---|
| 165 | l1norm[2] = l1_norm(edge[:,2]) |
|---|
| 166 | l2norm[0] = l2_norm(edge[:,0]) |
|---|
| 167 | l2norm[1] = l2_norm(edge[:,1]) |
|---|
| 168 | l2norm[2] = l2_norm(edge[:,2]) |
|---|
| 169 | linfnorm[0] = linf_norm(edge[:,0]) |
|---|
| 170 | linfnorm[1] = linf_norm(edge[:,1]) |
|---|
| 171 | linfnorm[2] = linf_norm(edge[:,2]) |
|---|
| 172 | l1list.append(l1norm) |
|---|
| 173 | l2list.append(l2norm) |
|---|
| 174 | linflist.append(linfnorm) |
|---|
| 175 | return (l1list, l2list, linflist) |
|---|
| 176 | |
|---|
| 177 | # Test an 8-way run of the shallow water equations |
|---|
| 178 | # against the sequential code. |
|---|
| 179 | |
|---|
| 180 | class Test_Parallel_Sw(unittest.TestCase): |
|---|
| 181 | def testParallelSw(self): |
|---|
| 182 | print "Expect this test to fail if not run from the parallel directory." |
|---|
| 183 | result = os.system("mpirun -np %d python test_parallel_sw.py" % nprocs) |
|---|
| 184 | assert_(result == 0) |
|---|
| 185 | |
|---|
| 186 | # Because we are doing assertions outside of the TestCase class |
|---|
| 187 | # the PyUnit defined assert_ function can't be used. |
|---|
| 188 | def assert_(condition, msg="Assertion Failed"): |
|---|
| 189 | if condition == False: |
|---|
| 190 | pypar.finalize() |
|---|
| 191 | raise AssertionError, msg |
|---|
| 192 | |
|---|
| 193 | if __name__=="__main__": |
|---|
| 194 | if pypar.size() == 1: |
|---|
| 195 | runner = unittest.TextTestRunner() |
|---|
| 196 | suite = unittest.makeSuite(Test_Parallel_Sw, 'test') |
|---|
| 197 | runner.run(suite) |
|---|
| 198 | else: |
|---|
| 199 | if pypar.rank() == 0: |
|---|
| 200 | l1norm_seq, l2norm_seq, linfnorm_seq = sequential_test() |
|---|
| 201 | l1norm_par, l2norm_par, linfnorm_par = parallel_test() |
|---|
| 202 | if pypar.rank() == 0: |
|---|
| 203 | assert_(len(l1norm_seq) == len(l1norm_par)) |
|---|
| 204 | assert_(len(l2norm_seq) == len(l2norm_par)) |
|---|
| 205 | assert_(len(linfnorm_seq) == len(linfnorm_par)) |
|---|
| 206 | assert_(len(l1norm_seq) == len(l2norm_seq)) |
|---|
| 207 | assert_(len(l2norm_seq) == len(linfnorm_seq)) |
|---|
| 208 | # Anything smaller than tol we consider to be 0. |
|---|
| 209 | tol = pow(10, -1 * double_precision()) |
|---|
| 210 | for x in range(len(l1norm_seq)): |
|---|
| 211 | for y in range(3): |
|---|
| 212 | assert_(abs(((l1norm_seq[x][y] - l1norm_par[x][y]) < tol))) |
|---|
| 213 | assert_(abs(((l2norm_seq[x][y] - l2norm_par[x][y]) < tol))) |
|---|
| 214 | assert_(abs(((linfnorm_seq[x][y] - linfnorm_par[x][y]) < tol))) |
|---|
| 215 | if x > 0: |
|---|
| 216 | assert_(abs(((l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol))) |
|---|
| 217 | assert_(abs(((l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol))) |
|---|
| 218 | assert_(abs(((linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol))) |
|---|
| 219 | assert_(abs(((l1norm_par[x][y] - l1norm_par[x-1][y]) < tol))) |
|---|
| 220 | assert_(abs(((l2norm_par[x][y] - l2norm_par[x-1][y]) < tol))) |
|---|
| 221 | assert_(abs(((linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol))) |
|---|
| 222 | pypar.finalize() |
|---|
| 223 | |
|---|
