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

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

Concentrating code

  • Property svn:executable set to *
File size: 8.3 KB
Line 
1
2"""Test a run of the sequential shallow water domain against
3a run of the parallel shallow water domain.
4
5WARNING: This assumes that the command to run jobs is mpirun.
6Tested with MPICH and LAM (Ole)
7"""
8
9#------------------------------------------------------------------------------
10# Import necessary modules
11#------------------------------------------------------------------------------
12
13import unittest
14import os
15import sys
16import pypar
17
18#from Numeric import allclose, array, zeros, Float, take, nonzero
19
20import numpy as num
21
22from anuga.pmesh.mesh_interface import create_mesh_from_regions
23
24from anuga.interface import rectangular_cross
25from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
26
27from anuga.utilities.numerical_tools import ensure_numeric
28from anuga.utilities.util_ext        import double_precision
29from anuga.utilities.norms           import l1_norm, l2_norm, linf_norm
30
31from anuga.interface import Domain
32from anuga.interface import Reflective_boundary
33from anuga.interface import Dirichlet_boundary
34from anuga.interface import Time_boundary
35from anuga.interface import Transmissive_boundary
36
37
38from anuga_parallel.parallel_api import distribute, myid, numprocs
39
40
41#--------------------------------------------------------------------------
42# Setup parameters
43#--------------------------------------------------------------------------
44
45mesh_filename = "merimbula_10785_1.tsh"
46#mesh_filename = "test-100.tsh"
47yieldstep = 1
48finaltime = 20
49quantity = 'stage'
50nprocs = 4
51
52#--------------------------------------------------------------------------
53# Setup procedures
54#--------------------------------------------------------------------------
55class 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
63
64    def __call__(self, x, y):
65        return self.h*((x>self.x0)&(x<self.x1))
66
67#--------------------------------------------------------------------------
68# Setup test
69#--------------------------------------------------------------------------
70def evolution_test(parallel=False):
71
72
73    domain = pmesh_to_domain_instance(mesh_filename, Domain)
74    domain.set_quantity('stage', Set_Stage(756000.0, 756500.0, 2.0))
75
76    #--------------------------------------------------------------------------
77    # Create parallel domain if requested
78    #--------------------------------------------------------------------------
79
80    if parallel:
81        if myid == 0: print 'DISTRIBUTING PARALLEL DOMAIN'
82        domain = distribute(domain, verbose=False)
83
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
90
91    domain.set_boundary({'outflow' :Br, 'inflow' :Br, 'inner' :Br, 'exterior' :Br, 'open' :Br})
92
93    #------------------------------------------------------------------------------
94    # Setup diagnostic arrays
95    #------------------------------------------------------------------------------
96    l1list = []
97    l2list = []
98    linflist = []
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)
103
104    #------------------------------------------------------------------------------
105    # Evolution
106    #------------------------------------------------------------------------------
107    if parallel:
108        if myid == 0: print 'PARALLEL EVOLVE'
109    else:
110        print 'SEQUENTIAL EVOLVE'
111       
112    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
113        edges = domain.quantities[quantity].edge_values.take(num.flatnonzero(domain.tri_full_flag),axis=0)
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])
120        linfnorm[0] = linf_norm(edges[:,0])
121        linfnorm[1] = linf_norm(edges[:,1])
122        linfnorm[2] = linf_norm(edges[:,2])
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:
128                domain.write_time()
129
130                #print edges[:,1]           
131                for p in range(1, numprocs):
132                    recv_norm = pypar.receive(p)
133                    l1norm += recv_norm
134                    recv_norm = pypar.receive(p)
135                    l2norm += recv_norm
136                    recv_norm = pypar.receive(p)
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:
153            domain.write_time()
154            l1list.append(l1norm)               
155            l2list.append(l2norm)
156            linflist.append(linfnorm)
157           
158
159    return (l1list, l2list, linflist)
160
161# Test an 8-way run of the shallow water equations
162# against the sequential code.
163
164class Test_Parallel_Sw(unittest.TestCase):
165    def testParallelSw(self):
166        print "Expect this test to fail if not run from the parallel directory."
167        result = os.system("mpirun -np %d python test_parallel_sw.py" % nprocs)
168        assert_(result == 0)
169
170# Because we are doing assertions outside of the TestCase class
171# the PyUnit defined assert_ function can't be used.
172def assert_(condition, msg="Assertion Failed"):
173    if condition == False:
174        #pypar.finalize()
175        raise AssertionError, msg
176
177if __name__=="__main__":
178    if numprocs == 1: 
179        runner = unittest.TextTestRunner()
180        suite = unittest.makeSuite(Test_Parallel_Sw, 'test')
181        runner.run(suite)
182    else:
183
184        pypar.barrier()
185        if myid == 0:
186            print 'SEQUENTIAL START'
187            l1norm_seq, l2norm_seq, linfnorm_seq = evolution_test(parallel=False)
188
189        pypar.barrier()
190        if myid ==0:
191            print 'PARALLEL START'
192       
193        l1norm_par, l2norm_par, linfnorm_par = evolution_test(parallel=True)
194       
195        if myid == 0:
196            assert_(len(l1norm_seq) == len(l1norm_par))
197            assert_(len(l2norm_seq) == len(l2norm_par))
198            assert_(len(linfnorm_seq) == len(linfnorm_par))
199            assert_(len(l1norm_seq) == len(l2norm_seq))
200            assert_(len(l2norm_seq) == len(linfnorm_seq))
201            # Anything smaller than tol we consider to be 0.
202            # This usualy comes out to be 10^-14 (DBL_DIG = 10^-15)
203            # * 10 to cater for rounding error in computation.
204            tol = pow(10, -1 * (double_precision() - 1))
205            for x in range(len(l1norm_seq)):
206                for y in range(3):
207                    # Calculate relative difference in the norms
208                    assert_(abs(l1norm_seq[x][y] - l1norm_par[x][y])/l1norm_seq[x][y] < tol)
209                    assert_(abs(l2norm_seq[x][y] - l2norm_par[x][y])/l2norm_seq[x][y] < tol)
210                    assert_(abs(linfnorm_seq[x][y] - linfnorm_par[x][y])/linfnorm_seq[x][y] < tol)
211                    if x > 0:
212                        # Verify that the quantity is being conserved across iterations.
213                        assert_(abs(l1norm_seq[x][y] - l1norm_seq[x-1][y]) < tol)
214                        assert_(abs(l2norm_seq[x][y] - l2norm_seq[x-1][y]) < tol)
215                        assert_(abs(linfnorm_seq[x][y] - linfnorm_seq[x-1][y]) < tol)
216                        assert_(abs(l1norm_par[x][y] - l1norm_par[x-1][y]) < tol)
217                        assert_(abs(l2norm_par[x][y] - l2norm_par[x-1][y]) < tol)
218                        assert_(abs(linfnorm_par[x][y] - linfnorm_par[x-1][y]) < tol)
219               
220            print 'Parallel test OK'
221
222
Note: See TracBrowser for help on using the repository browser.