source: trunk/anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py @ 7877

Last change on this file since 7877 was 7877, checked in by hudson, 14 years ago

Moved all development files into trunk.

  • Property svn:executable set to *
File size: 8.2 KB
RevLine 
[2738]1
[6047]2"""Test a run of the sequential shallow water domain against
3a run of the parallel shallow water domain.
[2738]4
[6047]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
[2738]13import unittest
14import os
15import sys
[6049]16import pypar
[2738]17
[7447]18import numpy as num
19
[6049]20
[2738]21
[7449]22
[6049]23from anuga.utilities.numerical_tools import ensure_numeric
24from anuga.utilities.util_ext        import double_precision
25from anuga.utilities.norms           import l1_norm, l2_norm, linf_norm
26
[7877]27from anuga import Domain
28from anuga import Reflective_boundary
29from anuga import Dirichlet_boundary
30from anuga import Time_boundary
31from anuga import Transmissive_boundary
[6047]32
[7877]33from anuga import rectangular_cross
34from anuga import create_domain_from_file
[6047]35
36
[7562]37from anuga_parallel.interface import distribute, myid, numprocs, finalize
[6047]38
[7449]39
[6047]40#--------------------------------------------------------------------------
41# Setup parameters
42#--------------------------------------------------------------------------
43
44mesh_filename = "merimbula_10785_1.tsh"
45#mesh_filename = "test-100.tsh"
46yieldstep = 1
47finaltime = 20
48quantity = 'stage'
[6048]49nprocs = 4
[7505]50verbose = False
[6047]51
52#--------------------------------------------------------------------------
53# Setup procedures
54#--------------------------------------------------------------------------
[2738]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
[6047]67#--------------------------------------------------------------------------
68# Setup test
69#--------------------------------------------------------------------------
70def 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]164class 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.
173def assert_(condition, msg="Assertion Failed"):
174    if condition == False:
[6721]175        #pypar.finalize()
[2776]176        raise AssertionError, msg
177
[2738]178if __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()
Note: See TracBrowser for help on using the repository browser.