source: development/parallel/test_parallel_sw.py @ 3460

Last change on this file since 3460 was 3460, checked in by jack, 18 years ago

Moved parallel to development from inundation

  • Property svn:executable set to *
File size: 8.9 KB
Line 
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"
8mesh_filename= "merimbula_10785_1.tsh"
9yieldstep = 1
10finaltime = 90
11quantity = 'stage'
12nprocs = 8
13
14import unittest
15import os
16import sys
17import print_stats
18import pypar
19
20from Numeric import array, zeros, Float, take, nonzero
21from pyvolution.shallow_water import Domain
22from pyvolution.shallow_water import Reflective_boundary as sw_reflective_boundary
23from pyvolution.shallow_water import Transmissive_boundary as sw_transmissive_boundary
24from parallel_shallow_water import Parallel_Domain
25from parallel_shallow_water import Reflective_boundary as par_reflective_boundary
26from parallel_shallow_water import Transmissive_boundary as par_transmissive_boundary
27from pyvolution.pmesh2domain import pmesh_to_domain_instance
28from utilities.norms import *
29from utilities.util_ext import double_precision
30from print_stats import print_test_stats, build_full_flag
31from pmesh_divide import pmesh_divide_metis
32from build_submesh import build_submesh
33from build_local import build_local_mesh
34from build_commun import send_submesh, rec_submesh, extract_hostmesh
35
36class 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
44
45    def __call__(self, x, y):
46        return self.h*((x>self.x0)&(x<self.x1))
47
48sys.path.append(".." + os.sep + "pyvolution")
49
50def 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
140def 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
179class 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.
187def assert_(condition, msg="Assertion Failed"):
188    if condition == False:
189        pypar.finalize()
190        raise AssertionError, msg
191
192if __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
Note: See TracBrowser for help on using the repository browser.