source: trunk/anuga_core/source/anuga_parallel/run_parallel_merimbula_test.py @ 8283

Last change on this file since 8283 was 3591, checked in by linda, 18 years ago

Fixed some bugs in build_submesh

File size: 4.5 KB
Line 
1#!/usr/bin/env python
2#########################################################
3#
4#  Main file for parallel mesh testing. Runs an advection
5# flow simulation using a rectangular mesh
6#
7#  This is a modification of the run_parallel_advection.py
8# file
9#
10#
11# *) The (new) files that have been added to manage the
12# grid partitioning are
13#    +) pmesh_divide_metis.py: subdivide a pmesh
14#    +) build_submesh.py: build the submeshes on the host
15# processor.
16#    +) build_local.py: build the GA mesh datastructure
17# on each processor.
18#    +) build_commun.py: handle the communication between
19# the host and processors
20#
21#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
22# June 2005
23#
24#
25#
26#########################################################
27
28import pypar    # The Python-MPI interface
29import time
30
31from sys import path
32
33print path
34
35# Numeric arrays
36from Numeric import array, zeros, Float
37
38from print_stats import print_test_stats, build_full_flag
39
40from anuga.abstract_2d_finite_volumes.pmesh2domain\
41     import pmesh_to_domain_instance
42from anuga.advection.advection import Domain as Advection_Domain
43from parallel_advection import Parallel_Domain
44
45from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
46     import Transmissive_boundary
47
48# mesh partition routines
49
50from pmesh_divide  import pmesh_divide_metis
51from build_submesh import build_submesh
52from build_local   import build_local_mesh
53from build_commun  import send_submesh, rec_submesh, extract_hostmesh
54
55
56
57class Set_Stage:
58    """Set an initial condition with constant water height, for x<x0
59    """
60
61    def __init__(self, x0=0.25, x1=0.5, h=1.0):
62        self.x0 = x0
63        self.x1 = x1
64        self.= h
65
66    def __call__(self, x, y):
67        return self.h*((x>self.x0)&(x<self.x1))
68
69# read in the processor information
70
71numprocs = pypar.size()
72myid = pypar.rank()
73processor_name = pypar.Get_processor_name()
74
75
76#-------
77# Domain
78rect = zeros( 4, Float) # Buffer for results
79
80if myid == 0:
81
82    # read in the test files
83
84    filename = 'test-100.tsh'
85#    filename = 'merimbula_10785.tsh'
86    nx = numprocs
87    ny = 1
88    if nx*ny != numprocs:
89        print "WARNING: number of subboxes is not equal to the number of proc"
90
91    domain_full = pmesh_to_domain_instance(filename, Advection_Domain)
92    domain_full.set_quantity('stage', Set_Stage(200.0,300.0,1.0))
93#    domain_full.set_quantity('stage', Set_Stage(756000.0,756500.0,4.0))
94
95    nodes, triangles, boundary, triangles_per_proc, quantities  =\
96            pmesh_divide_metis(domain_full, numprocs)
97   
98    # subdivide the mesh
99    rect = array(domain_full.xy_extent, Float)
100#    rect = array(rect, Float)
101
102    submesh = build_submesh(nodes, triangles, boundary, quantities, \
103                            triangles_per_proc)
104
105    # send the mesh partition to the appropriate processor
106
107    for p in range(1, numprocs):
108      send_submesh(submesh, triangles_per_proc, p)
109
110    # Build the local mesh for processor 0
111   
112    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
113              extract_hostmesh(submesh, triangles_per_proc)
114
115# read in the mesh partition that belongs to this
116# processor (note that the information is in the
117# correct form for the GA data structure
118
119else:
120    [points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict] = \
121             rec_submesh(0)
122
123pypar.broadcast(rect,0)
124
125domain = Parallel_Domain(points, vertices, boundary,
126                                   full_send_dict  = full_send_dict,
127                                   ghost_recv_dict = ghost_recv_dict,
128                                   velocity = [0.1,0.0])
129
130# Make a notes of which triangles are full and which are ghost
131
132tri_full_flag = build_full_flag(domain, ghost_recv_dict)
133
134#domain.initialise_visualiser(rect=rect)
135
136#Boundaries
137
138T = Transmissive_boundary(domain)
139#R = Reflective_boundary(domain)
140domain.set_boundary( {'outflow': T, 'inflow': T, 'inner':T, 'exterior': T, 'open':T, 'ghost':None} )
141
142
143
144domain.set_quantity('stage', quantities['stage'])
145
146#---------
147# Evolution
148t0 = time.time()
149
150
151try:
152    domain.initialise_visualiser(rect=rect)
153    #domain.visualiser.coloring['stage'] = True
154    #domain.visualiser.scale_z['stage'] = 0.2
155except:
156    print 'No visualiser'
157
158
159from norms import linf_norm
160
161yieldstep = 1
162finaltime = 200
163
164#yieldstep = 1000
165#finaltime = 50000
166
167for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
168    if myid == 0:
169        domain.write_time()
170    print_test_stats(domain, tri_full_flag)
171
172if myid == 0:
173    print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.