source: anuga_core/source/anuga_parallel/test_parallel_sw.py @ 6048

Last change on this file since 6048 was 6048, checked in by steve, 15 years ago

Just madeteh default number of processors to be 4 in test_parallel_sw.py

  • Property svn:executable set to *
File size: 8.4 KB
Line 
1#!/usr/bin/env python
2
3"""Test a run of the sequential shallow water domain against
4a run of the parallel shallow water domain.
5
6WARNING: This assumes that the command to run jobs is mpirun.
7Tested with MPICH and LAM (Ole)
8"""
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14import unittest
15import os
16import sys
17import print_stats
18
19from Numeric import allclose, array, zeros, Float, take, nonzero
20
21from anuga.pmesh.mesh_interface import create_mesh_from_regions
22from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
23from anuga.utilities.numerical_tools import ensure_numeric
24from anuga.utilities.polygon import is_inside_polygon
25from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
26from anuga.utilities.util_ext import double_precision
27from anuga.utilities.norms import l1_norm, l2_norm, linf_norm
28
29from anuga.shallow_water import Domain
30from anuga.shallow_water import Reflective_boundary
31from anuga.shallow_water import Dirichlet_boundary
32from anuga.shallow_water import Time_boundary
33from anuga.shallow_water import Transmissive_boundary
34
35import pypar
36
37from parallel_api import distribute, myid, numprocs
38
39
40#--------------------------------------------------------------------------
41# Setup parameters
42#--------------------------------------------------------------------------
43
44mesh_filename = "merimbula_10785_1.tsh"
45#mesh_filename = "test-100.tsh"
46yieldstep = 1
47finaltime = 20
48quantity = 'stage'
49nprocs = 4
50
51#--------------------------------------------------------------------------
52# Setup procedures
53#--------------------------------------------------------------------------
54class 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
62
63    def __call__(self, x, y):
64        return self.h*((x>self.x0)&(x<self.x1))
65
66#--------------------------------------------------------------------------
67# Setup test
68#--------------------------------------------------------------------------
69def 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
162class 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        print 'result ',result
167        assert_(result == 0)
168
169# Because we are doing assertions outside of the TestCase class
170# the PyUnit defined assert_ function can't be used.
171def assert_(condition, msg="Assertion Failed"):
172    if condition == False:
173        pypar.finalize()
174        raise AssertionError, msg
175
176if __name__=="__main__":
177    if numprocs == 1: 
178        runner = unittest.TextTestRunner()
179        suite = unittest.makeSuite(Test_Parallel_Sw, 'test')
180        runner.run(suite)
181    else:
182
183        pypar.barrier()
184        if myid == 0:
185            print 'SEQUENTIAL START'
186            l1norm_seq, l2norm_seq, linfnorm_seq = evolution_test(parallel=False)
187
188        pypar.barrier()
189        if myid ==0:
190            print 'PARALLEL START'
191       
192        l1norm_par, l2norm_par, linfnorm_par = evolution_test(parallel=True)
193       
194        if myid == 0:
195            assert_(len(l1norm_seq) == len(l1norm_par))
196            assert_(len(l2norm_seq) == len(l2norm_par))
197            assert_(len(linfnorm_seq) == len(linfnorm_par))
198            assert_(len(l1norm_seq) == len(l2norm_seq))
199            assert_(len(l2norm_seq) == len(linfnorm_seq))
200            # Anything smaller than tol we consider to be 0.
201            # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15)
202            # * 10 to cater for rounding error in computation.
203            tol = pow(10, -1 * (double_precision() - 1))
204            for x in range(len(l1norm_seq)):
205                for y in range(3):
206                    # Calculate relative difference in the norms
207                    assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol)
208                    assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol)
209                    assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol)
210                    if x > 0:
211                        # Verify that the quantity is being conserved across iterations.
212                        assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol)
213                        assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol)
214                        assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol)
215                        assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol)
216                        assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol)
217                        assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol)
218
219        if myid == 0:
220            print 'Parallel test OK'
221
222
Note: See TracBrowser for help on using the repository browser.