source: anuga_work/development/parallel/test_parallel_sw.py @ 5175

Last change on this file since 5175 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

  • 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 anuga.pyvolution.shallow_water import Domain
22from anuga.pyvolution.shallow_water import Reflective_boundary as sw_reflective_boundary
23from anuga.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 anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
28from anuga.utilities.norms import *
29from anuga.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.