source: trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py @ 8262

Last change on this file since 8262 was 8262, checked in by steve, 12 years ago

Added in unit test for dump_triangulation

File size: 6.6 KB
Line 
1"""Class Parallel_shallow_water_domain -
22D triangular domains for finite-volume computations of
3the shallow water equation, with extra structures to allow
4communication between other Parallel_domains and itself
5
6This module contains a specialisation of class Domain
7from module shallow_water.py
8
9Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
10Geoscience Australia, 2004-2005
11
12"""
13
14from anuga import Domain
15
16from anuga_parallel.parallel_generic_communications import *
17from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
18
19import numpy as num
20
21
22#Import matplotlib
23
24
25
26class Parallel_domain(Domain):
27
28    def __init__(self, coordinates, vertices,
29                 boundary=None,
30                 full_send_dict=None,
31                 ghost_recv_dict=None,
32                 number_of_full_nodes=None,
33                 number_of_full_triangles=None,
34                 geo_reference=None,
35                 tri_map=None,
36                 inv_tri_map=None): #jj added this
37
38        Domain.__init__(self,
39                        coordinates,
40                        vertices,
41                        boundary,
42                        full_send_dict=full_send_dict,
43                        ghost_recv_dict=ghost_recv_dict,
44                        processor=pypar.rank(),
45                        numproc=pypar.size(),
46                        number_of_full_nodes=number_of_full_nodes,
47                        number_of_full_triangles=number_of_full_triangles,
48                        geo_reference=geo_reference) #jj added this
49       
50        # PETE: Find the number of full nodes and full triangles, this is a temporary fix
51        # until the bug with get_number_of_full_[nodes|triangles]() is fixed.
52
53        if number_of_full_nodes is not None:
54            self.number_of_full_nodes_tmp = number_of_full_nodes
55        else:
56            self.number_of_full_nodes_tmp = get_number_of_nodes()
57
58        if number_of_full_triangles is not None:
59            self.number_of_full_triangles_tmp = number_of_full_triangles
60        else:
61            self.number_of_full_triangles_tmp = get_number_of_triangles()
62
63        setup_buffers(self)
64
65        self.tri_map = tri_map
66        self.inv_tri_map = inv_tri_map
67
68
69    def set_name(self, name):
70        """Assign name based on processor number
71        """
72
73        if name.endswith('.sww'):
74            name = name[:-4]
75
76        # Call parents method with processor number attached.
77        Domain.set_name(self, name + '_P%d_%d' %(self.processor, self.numproc))
78
79
80    def update_timestep(self, yieldstep, finaltime):
81        """Calculate local timestep
82        """
83
84        communicate_flux_timestep(self, yieldstep, finaltime)
85
86        Domain.update_timestep(self, yieldstep, finaltime)
87
88
89
90    def update_ghosts(self):
91
92        # We must send the information from the full cells and
93        # receive the information for the ghost cells
94
95
96        communicate_ghosts(self)
97
98    def apply_fractional_steps(self):
99
100        for operator in self.fractional_step_operators:
101            operator()
102
103        # PETE: Make sure that there are no deadlocks here
104
105        self.update_ghosts()
106
107# =======================================================================
108# PETE: NEW METHODS FOR FOR PARALLEL STRUCTURES. Note that we assume the
109# first "number_of_full_[nodes|triangles]" are full [nodes|triangles]
110# For full triangles it is possible to enquire self.tri_full_flag == True
111# =======================================================================
112
113    def get_number_of_full_triangles(self, *args, **kwargs):
114        return self.number_of_full_triangles_tmp
115
116    def get_full_centroid_coordinates(self, *args, **kwargs):
117        C = self.mesh.get_centroid_coordinates(*args, **kwargs)
118        return C[:self.number_of_full_triangles_tmp, :]
119
120    def get_full_vertex_coordinates(self, *args, **kwargs):
121        V = self.mesh.get_vertex_coordinates(*args, **kwargs)
122        return V[:3*self.number_of_full_triangles_tmp,:]
123
124    def get_full_triangles(self, *args, **kwargs):
125        T = self.mesh.get_triangles(*args, **kwargs)
126        return T[:self.number_of_full_triangles_tmp,:]
127
128    def get_full_nodes(self, *args, **kwargs):
129        N = self.mesh.get_nodes(*args, **kwargs)
130        return N[:self.number_of_full_nodes_tmp,:]
131
132    def get_tri_map(self):
133        return self.tri_map
134
135    def get_inv_tri_map(self):
136        return self.inv_tri_map
137
138    '''
139    Outputs domain triangulation, full triangles are shown in green while ghost triangles are shown in blue.
140    The default filename is "domain.png"
141    '''
142    def dump_triangulation(self, filename="domain.png"):
143        # Get vertex coordinates, partition full and ghost triangles based on self.tri_full_flag
144
145        try:
146            import matplotlib
147            matplotlib.use('Agg')
148            import matplotlib.pyplot as plt
149            import matplotlib.tri as tri
150        except:
151            print "Couldn't import module from matplotlib, probably you need to update matplotlib"
152            raise
153
154        vertices = self.get_vertex_coordinates()
155        full_mask = num.repeat(self.tri_full_flag == 1, 3)
156        ghost_mask = num.repeat(self.tri_full_flag == 0, 3)
157       
158        myid = pypar.rank()
159        numprocs = pypar.size()
160
161        if myid == 0:
162            fx = {}
163            fy = {}
164            gx = {}
165            gy = {}
166
167            # Proc 0 gathers full and ghost nodes from self and other processors
168            fx[0] = vertices[full_mask,0]
169            fy[0] = vertices[full_mask,1]
170            gx[0] = vertices[ghost_mask,0]
171            gy[0] = vertices[ghost_mask,1]
172           
173            for i in range(1,numprocs):
174                fx[i] = pypar.receive(i)
175                fy[i] = pypar.receive(i)
176                gx[i] = pypar.receive(i)
177                gy[i] = pypar.receive(i)
178
179            # Plot full triangles
180            for i in range(0, numprocs):
181                n = int(len(fx[i])/3)
182                           
183                triang = num.array(range(0,3*n))
184                triang.shape = (n, 3)
185                plt.triplot(fx[i], fy[i], triang, 'g-')
186
187            # Plot ghost triangles
188            for i in range(0, numprocs):
189                n = int(len(gx[i])/3)
190                           
191                triang = num.array(range(0,3*n))
192                triang.shape = (n, 3)
193                plt.triplot(gx[i], gy[i], triang, 'b--')
194
195            # Save triangulation to location pointed by filename
196            plt.savefig(filename)
197
198        else:
199            # Proc 1..numprocs send full and ghost triangles to Proc 0
200            pypar.send(vertices[full_mask,0], 0)
201            pypar.send(vertices[full_mask,1], 0)
202            pypar.send(vertices[ghost_mask,0], 0)
203            pypar.send(vertices[ghost_mask,1], 0)
204
Note: See TracBrowser for help on using the repository browser.