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