# Changeset 7400

Ignore:
Timestamp:
Aug 21, 2009, 2:02:46 PM (15 years ago)
Message:

Commit a working copy of numpy version of build_commun

Location:
anuga_core/source/anuga_parallel
Files:
11 edited

Unmodified
Removed

• ## anuga_core/source/anuga_parallel/build_local.py

 r3579 #  Authors: Linda Stals and Matthew Hardy, June 2005 #  Modified: Linda Stals, Nov 2005 (optimise python code) # # ######################################################### from Numeric import  zeros, Float, Int, concatenate, \ take, arrayrange, put, sort, compress, equal #            Steve Roberts, Aug 2009 (updating to numpy) # # ######################################################### #from Numeric import  zeros, Float, Int, concatenate, \ #     take, arrayrange, put, sort, compress, equal import numpy as num ######################################################### # Convert the format of the data to that used by ANUGA # Extract the nodes (using the local ID) GAnodes = take(nodes, (1, 2), 1) GAnodes = num.take(nodes, (1, 2), 1) # Build a global ID to local ID mapping if nodes[i][0] > NGlobal: NGlobal = nodes[i][0] index = zeros(int(NGlobal)+1, Int) put(index, take(nodes, (0,), 1).astype(Int), \ arrayrange(Nnodes)) index = num.zeros(int(NGlobal)+1, num.int) num.put(index, num.take(nodes, (0,), 1).astype(num.int), \ num.arange(Nnodes)) # Change the global IDs in the triangles to the local IDs GAtriangles = zeros((Ntriangles, 3), Int) GAtriangles[:,0] = take(index, triangles[:,0]) GAtriangles[:,1] = take(index, triangles[:,1]) GAtriangles[:,2] = take(index, triangles[:,2]) GAtriangles = num.zeros((Ntriangles, 3), num.int) GAtriangles[:,0] = num.take(index, triangles[:,0]) GAtriangles[:,1] = num.take(index, triangles[:,1]) GAtriangles[:,2] = num.take(index, triangles[:,2]) # Change the triangle numbering in the boundaries # information by the global numbering) ghostc = sort(ghostc, 0) ghostc = num.sort(ghostc, 0) for c in range(nproc): s = ghostc[:,0] d = compress(equal(ghostc[:,1],c), s) d = num.compress(num.equal(ghostc[:,1],c), s) if len(d) > 0: ghost_recv[c] = [0, 0] ghost_recv[c][1] = d ghost_recv[c][0] = take(index, d) ghost_recv[c][0] = num.take(index, d) # Build a temporary copy of the full_send dictionary for neigh in tmp_send: neigh_commun = sort(tmp_send[neigh], 0) neigh_commun = num.sort(tmp_send[neigh], 0) full_send[neigh] = [0, 0] full_send[neigh][0] = neigh_commun[:,1] # Combine the full nodes and ghost nodes nodes = concatenate((submesh["full_nodes"], \ nodes = num.concatenate((submesh["full_nodes"], \ submesh["ghost_nodes"])) # Combine the full triangles and ghost triangles gtri =  take(submesh["ghost_triangles"],(1, 2, 3),1) triangles = concatenate((submesh["full_triangles"], gtri)) gtri =  num.take(submesh["ghost_triangles"],(1, 2, 3),1) triangles = num.concatenate((submesh["full_triangles"], gtri)) # Combine the full boundaries and ghost boundaries if id > NGlobal: NGlobal = id index = zeros(int(NGlobal)+1, Int) index[lower_t:upper_t]=arrayrange(upper_t-lower_t) index = num.zeros(int(NGlobal)+1, num.int) index[lower_t:upper_t]=num.arange(upper_t-lower_t) for i in range(len(submesh["ghost_triangles"])): index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t Nf = len(submesh["full_quan"][k]) Ng = len(submesh["ghost_quan"][k]) quantities[k] = zeros((Nf+Ng, 3), Float) quantities[k] = num.zeros((Nf+Ng, 3), num.float) quantities[k][0:Nf] = submesh["full_quan"][k] quantities[k][Nf:Nf+Ng] = submesh["ghost_quan"][k] return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \ full_send
• ## anuga_core/source/anuga_parallel/build_submesh.py

 r7389 ######################################################### # # Subdivide the GA domain. This module is primarily # Subdivide the domain. This module is primarily # responsible for building the ghost layer and # communication pattern #  Author: Linda Stals, June 2005 #  Modified: Linda Stals, Nov 2005 (optimise python code) #            Steve Roberts, Aug 2009 (convert to numpy) # # import sys from Numeric import zeros, Float, Int, concatenate, \ reshape, arrayrange, take, nonzero #from Numeric import zeros, Float, Int, concatenate, \ #     reshape, arrayrange, take, nonzero import numpy as num from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh boundary_list = [] submesh = {} node_range = reshape(arrayrange(nnodes),(nnodes,1)) tsubnodes = concatenate((node_range, nodes), 1) node_range = num.reshape(num.arange(nnodes),(nnodes,1)) #print node_range tsubnodes = num.concatenate((node_range, nodes), 1) # Loop over processors # Find nodes in processor p nodemap = zeros(nnodes, 'i') nodemap = num.zeros(nnodes, 'i') for t in subtriangles: nodemap[t[0]]=1 nodemap[t[2]]=1 node_list.append(take(tsubnodes,nonzero(nodemap))) node_list.append(tsubnodes.take(num.flatnonzero(nodemap),axis=0)) # Move to the next processor # Find the first layer of boundary triangles trianglemap = zeros(ntriangles, 'i') trianglemap = num.zeros(ntriangles, 'i') for t in range(tlower, tupper): n = mesh.neighbours[t, 0] if n >= 0: if n < tlower or n >= tupper: # Build the triangle list and make note of the vertices nodemap = zeros(ncoord, 'i') nodemap = num.zeros(ncoord, 'i') fullnodes = submesh["full_nodes"][p] nodemap[t[2]] = 1 trilist = reshape(arrayrange(ntriangles),(ntriangles,1)) tsubtriangles = concatenate((trilist, mesh.triangles), 1) subtriangles = take(tsubtriangles, nonzero(trianglemap)) trilist = num.reshape(num.arange(ntriangles),(ntriangles,1)) tsubtriangles = num.concatenate((trilist, mesh.triangles), 1) subtriangles = tsubtriangles.take(num.flatnonzero(trianglemap),axis=0) # Keep a record of the triangle vertices, if they are not already there nodemap[int(n[0])] = 0 nodelist = reshape(arrayrange(ncoord),(ncoord,1)) tsubnodes = concatenate((nodelist, mesh.get_nodes()), 1) subnodes = take(tsubnodes, nonzero(nodemap)) nodelist = num.reshape(num.arange(ncoord),(ncoord,1)) tsubnodes = num.concatenate((nodelist, mesh.get_nodes()), 1) subnodes = tsubnodes.take(num.flatnonzero(nodemap),axis=0) # Clean up before exiting ######################################################### def is_in_processor(ghost_list, tlower, tupper, n): return (n in ghost_list) or (tlower <= n and tupper > n) return num.equal(ghost_list,n).any() or (tlower <= n and tupper > n) def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p): ghost_list = [] subboundary = {} for t in ghosttri: ghost_list.append(t[0]) for t in ghosttri: n = mesh.neighbours[t[0], 0] if not is_in_processor(ghost_list, tlower, tupper, n): # Loop over the ghost triangles ghost_commun = zeros((len(subtri), 2), Int) ghost_commun = num.zeros((len(subtri), 2), num.int) for i in range(len(subtri)): ghost_nodes.append(subnodes) # Find the boundary layer formed by the ghost triangles for k in quantities: submesh["full_quan"][k].append(quantities[k][lower:upper]) submesh["ghost_quan"][k].append(zeros( (M,3) , Float)) submesh["ghost_quan"][k].append(num.zeros( (M,3) , num.float)) for j in range(M): submesh["ghost_quan"][k][p][j] = \

 r5763 from anuga.advection import * from Numeric import zeros, Float, Int, ones, allclose, array #from Numeric import zeros, Float, Int, ones, allclose, array import numpy as num import pypar # For some reason it looks like pypar only reduces numeric arrays # hence we need to create some dummy arrays for communication ltimestep = ones( 1, Float ) ltimestep = num.ones( 1, num.float ) ltimestep[0] = self.flux_timestep gtimestep = zeros( 1, Float) # Buffer for results pypar.raw_reduce(ltimestep, gtimestep, pypar.MIN, 0) gtimestep = num.zeros( 1, num.float ) # Buffer for results #ltimestep = self.flux_timeste #print self.processor, ltimestep, gtimestep pypar.reduce(ltimestep, pypar.MIN, 0, buffer=gtimestep) #print self.processor, ltimestep, gtimestep pypar.broadcast(gtimestep,0) #print self.processor, ltimestep, gtimestep self.flux_timestep = gtimestep[0] # the separate processors from Numeric import take,put #from Numeric import take,put import numpy as num import time t0 = time.time() #for i in range(N): #    Xout[i,0] = stage_cv[Idf[i]] Xout[:,0] = take(stage_cv, Idf) Xout[:,0] = num.take(stage_cv, Idf) pypar.send(Xout,send_proc) N = len(Idg) put(stage_cv, Idg, X[:,0]) num.put(stage_cv, Idg, X[:,0]) #for i in range(N): #    stage_cv[Idg[i]] = X[i,0] #    stage_cv[Idg[i]] = stage_cv[Idf[i]] put(stage_cv, Idg, take(stage_cv, Idf)) num.put(stage_cv, Idg, num.take(stage_cv, Idf))
• ## anuga_core/source/anuga_parallel/parallel_api.py

 r6721 """ from Numeric import zeros # The abstract Python-MPI interface # Build the mesh that should be assigned to each processor, # this includes ghost nodes and the communication pattern
• ## anuga_core/source/anuga_parallel/parallel_meshes.py

 r6721 import sys from Numeric import array, zeros, Float, Int, ones, sum #from Numeric import array, zeros, Float, Int, ones, sum import numpy as num import pypar E = EIndex(n,m) points = zeros( (Np,2), Float) points = num.zeros( (Np,2), num.float) for i in range(m+1): elements = zeros( (Nt,3), Int) elements = num.zeros( (Nt,3), num.int) boundary = {} Idgl = [] print Idgr Idfl = array(Idfl,Int) Idgr = array(Idgr,Int) Idfl = num.array(Idfl,num.int) Idgr = num.array(Idgr,num.int) print Idfl Idfl.extend(Idfr) Idgr.extend(Idgl) Idfl = array(Idfl,Int) Idgr = array(Idgr,Int) Idfl = num.array(Idfl,num.int) Idgr = num.array(Idgr,num.int) full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl] ghost_recv_dict[(processor-1)%numproc] = [Idgr, Idgr] else: Idfl = array(Idfl,Int) Idgl = array(Idgl,Int) Idfr = array(Idfr,Int) Idgr = array(Idgr,Int) Idfl = num.array(Idfl,num.int) Idgl = num.array(Idgl,num.int) Idfr = num.array(Idfr,num.int) Idgr = num.array(Idgr,num.int) full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl] E = EIndex(n,m) points = zeros( (Np,2), Float) points = num.zeros( (Np,2), num.float) for i in range(m+1): elements = zeros( (Nt,3), Int) elements = num.zeros( (Nt,3), num.int) boundary = {} ghosts = {} E = EIndex(n,m) points = zeros( (Np,2), Float) points = num.zeros( (Np,2), num.float) for i in range(m+1): elements = zeros( (Nt,3), Int) elements = num.zeros( (Nt,3), num.int) boundary = {} ghosts = {}
• ## anuga_core/source/anuga_parallel/parallel_shallow_water.py

 r5763 from anuga.shallow_water.shallow_water_domain import * from Numeric import zeros, Float, Int, ones, allclose, array import numpy as num import pypar #        for key in full_send_dict: #            buffer_shape = full_send_dict[key][0].shape[0] #            full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) #            full_send_dict[key].append(num.zeros( (buffer_shape,self.nsys) ,num.loat)) # # #        for key in ghost_recv_dict: #            buffer_shape = ghost_recv_dict[key][0].shape[0] #            ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) #            ghost_recv_dict[key].append(num.zeros( (buffer_shape,self.nsys) ,num.float)) # #        self.full_send_dict  = full_send_dict # Buffers for synchronisation of timesteps self.local_timestep = zeros(1, Float) self.global_timestep = zeros(1, Float) self.local_timesteps = zeros(self.numproc, Float) self.local_timestep = num.zeros(1, num.float) self.global_timestep = num.zeros(1, num.float) self.local_timesteps = num.zeros(self.numproc, num.loat)
• ## anuga_core/source/anuga_parallel/pmesh_divide.py

 r6721 #  Modified: Linda Stals, Nov 2005 #            Jack Kelly, Nov 2005 #            Steve Roberts, Aug 2009 (updating to numpy) # # from math import floor from Numeric import zeros, Float, Int, reshape, argsort, ArrayType import numpy as num #import zeros, float, Int, reshape, argsort, ArrayType ######################################################### # Temporary storage area index = zeros(N, Int) index = num.zeros(N, num.int) q_reord = {} for k in quantities: q_reord[k] = zeros((N, 3), Float) q_reord[k] = num.zeros((N, 3), num.float) for i in range(N): q_reord[k][index[i]]=quantities[k].vertex_values[i] try: from pymetis.metis import partMeshNodal except ImportError: print "***************************************************" raise ImportError def pmesh_divide_metis(domain, n_procs): n_vert = domain.get_number_of_nodes() t_list = domain.triangles.copy() t_list = reshape(t_list, (-1,)) t_list = num.reshape(t_list, (-1,)) # The 1 here is for triangular mesh elements. # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero # dimensional arrays. Correct this. if type(epart[0]) == ArrayType: epart_new = zeros(len(epart), Int) if type(epart[0]) == num.ndarray: epart_new = num.zeros(len(epart), num.int) for i in range(len(epart)): epart_new[i] = epart[i][0] quantities = {} for k in domain.quantities: quantities[k] = zeros((n_tri, 3), Float) quantities[k] = num.zeros((n_tri, 3), num.float) for i in range(n_tri): quantities[k][i] = domain.quantities[k].vertex_values[i] # this helps with the communication ttriangles = zeros((len(triangles), 3), Int) ttriangles = num.zeros((len(triangles), 3), num.int) for i in range(len(triangles)): ttriangles[i] = triangles[i]