Changeset 1414


Ignore:
Timestamp:
May 17, 2005, 6:28:13 PM (20 years ago)
Author:
steve
Message:
 
Location:
inundation/ga/storm_surge
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/Hobart/run_hobart_buildings.py

    r1404 r1414  
    4242domain.default_order = 1
    4343domain.smooth = True
    44 domain.visualise = False
     44domain.visualise = True
    4545
    4646#------------------------------
  • inundation/ga/storm_surge/parallel/parallel_advection.py

    r1407 r1414  
    1616from advection import *
    1717Advection_Domain = Domain
    18 from Numeric import zeros, Float, Int, ones
     18from Numeric import zeros, Float, Int, ones, allclose,  array
    1919import pypar
    2020
     
    2222class Parallel_Domain(Advection_Domain):
    2323
    24     def __init__(self, coordinates, vertices, boundary = None, ghosts = None, velocity = None, processor=0):
     24    def __init__(self, coordinates, vertices, boundary = None, ghosts = None,
     25                 velocity = None):
    2526
    2627        Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
     
    2829        N = self.number_of_elements
    2930
    30         self.processor = processor
     31        self.processor = pypar.rank()
     32        self.numproc   = pypar.size()
    3133
    3234        if  ghosts == None:
     
    5254
    5355    def update_timestep(self, yieldstep, finaltime):
     56
     57        # Calculate local timestep
    5458        Advection_Domain.update_timestep(self, yieldstep, finaltime)
    5559
    56         timestep = 0.0
    57         pypar.raw_reduce(self.timestep,timestep,pypar.MIN,0)
    58         print 'Processor %d, Timestep %.4f'%(self.processor,timestep)
    59 
     60        # For some reason it looks like pypar only reduces numeric arrays
     61        # hence we need to create some dummy arrays for communication
     62        ltimestep = ones( 1, Float )
     63        ltimestep[0] = self.timestep
     64
     65        gtimestep = zeros( 1, Float) # Buffer for results
     66
     67        pypar.raw_reduce(ltimestep, gtimestep, pypar.MIN, 0)
     68        pypar.broadcast(gtimestep,0)
     69        pypar.Barrier()
     70
     71        self.timestep = gtimestep[0]
    6072
    6173
    6274
    6375    def update_ghosts(self):
     76
     77        # We must send the information from the full cells and
     78        # receive the information for the ghost cells
     79        # We have a dictionary of lists with ghosts expecting updates from
     80        # the separate processors
     81
     82        stage_cv = self.quantities['stage'].centroid_values
     83
     84#        for iproc in full_send_dict:
     85#            Ids = full_send_dict[iproc][0]
     86#            X   = full_send_dict[iproc][1]
     87#            for i _ in enumerate(X):
     88#                X[i] = stage_cv[Ids[i]]
     89#            pypar.send(X,iproc)
     90#
     91#        for iproc in ghost_recv_dict:
     92#            Ids = ghost_recv_dict[iproc][0]
     93#            X   = ghost_recv_dict[iproc][1]
     94#            X = pypar.receive(iproc,X)
     95#            for i _ in enumerate(X):
     96#                stage_cv[Ids[i]] = X[i]
     97
     98
    6499        if self.ghosts is not None:
    65100            stage_cv = self.quantities['stage'].centroid_values
     
    101136
    102137
    103 def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
     138
     139def Parallel_rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
    104140
    105141
     
    107143    with m+1 by n+1 grid points
    108144    and side lengths len1, len2. If side lengths are omitted
    109     the mesh defaults to the unit square.
     145    the mesh defaults to the unit square, divided between all the
     146    processors
    110147
    111148    len1: x direction (left to right)
    112149    len2: y direction (bottom to top)
    113150
    114     Return to lists: points and elements suitable for creating a Mesh or
    115     FVMesh object, e.g. Mesh(points, elements)
    116151    """
    117152
     
    178213                ghosts[nt] = E(m-2,j)
    179214
    180             if j == n-1:
    181                 ghosts[nt] = E(i,1)
    182 
    183             if j == 0:
    184                 ghosts[nt] = E(i,n-2)
    185 
    186215            if i == m-1:
    187216                boundary[nt, 2] = 'right'
     
    197226                ghosts[nt] = E(m-2,j)+1
    198227
    199             if j == n-1:
    200                 ghosts[nt] = E(i,1)+1
    201 
    202             if j == 0:
    203                 ghosts[nt] = E(i,n-2)+1
    204 
    205228            if i == 0:
    206229                boundary[nt, 2] = 'left'
     
    209232            elements[nt,:] = [i1,i2,i3]
    210233
    211     #bottom left
    212     nt = E(0,0)
    213     nf = E(m-2,n-2)
    214     ghosts[nt]   = nf
    215     ghosts[nt+1] = nf+1
    216 
    217     #bottom right
    218     nt = E(m-1,0)
    219     nf = E(1,n-2)
    220     ghosts[nt]   = nf
    221     ghosts[nt+1] = nf+1
    222 
    223     #top left
    224     nt = E(0,n-1)
    225     nf = E(m-2,1)
    226     ghosts[nt]   = nf
    227     ghosts[nt+1] = nf+1
    228 
    229     #top right
    230     nt = E(m-1,n-1)
    231     nf = E(1,1)
    232     ghosts[nt]   = nf
    233     ghosts[nt+1] = nf+1
    234234
    235235    return points, elements, boundary, ghosts
    236236
    237 def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
     237def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
    238238
    239239
     
    312312                ghosts[nt] = E(m-2,j)
    313313
     314            if j == n-1:
     315                ghosts[nt] = E(i,1)
     316
     317            if j == 0:
     318                ghosts[nt] = E(i,n-2)
     319
    314320            if i == m-1:
    315321                boundary[nt, 2] = 'right'
     
    325331                ghosts[nt] = E(m-2,j)+1
    326332
     333            if j == n-1:
     334                ghosts[nt] = E(i,1)+1
     335
     336            if j == 0:
     337                ghosts[nt] = E(i,n-2)+1
     338
    327339            if i == 0:
    328340                boundary[nt, 2] = 'left'
     
    331343            elements[nt,:] = [i1,i2,i3]
    332344
     345    #bottom left
     346    nt = E(0,0)
     347    nf = E(m-2,n-2)
     348    ghosts[nt]   = nf
     349    ghosts[nt+1] = nf+1
     350
     351    #bottom right
     352    nt = E(m-1,0)
     353    nf = E(1,n-2)
     354    ghosts[nt]   = nf
     355    ghosts[nt+1] = nf+1
     356
     357    #top left
     358    nt = E(0,n-1)
     359    nf = E(m-2,1)
     360    ghosts[nt]   = nf
     361    ghosts[nt+1] = nf+1
     362
     363    #top right
     364    nt = E(m-1,n-1)
     365    nf = E(1,1)
     366    ghosts[nt]   = nf
     367    ghosts[nt+1] = nf+1
    333368
    334369    return points, elements, boundary, ghosts
     370
     371def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
     372
     373
     374    """Setup a rectangular grid of triangles
     375    with m+1 by n+1 grid points
     376    and side lengths len1, len2. If side lengths are omitted
     377    the mesh defaults to the unit square.
     378
     379    len1: x direction (left to right)
     380    len2: y direction (bottom to top)
     381
     382    Return to lists: points and elements suitable for creating a Mesh or
     383    Domain object, e.g. Mesh(points, elements)
     384    """
     385
     386    from config import epsilon
     387    from Numeric import zeros, Float, Int
     388
     389    delta1 = float(len1)/m
     390    delta2 = float(len2)/n
     391
     392    #Calculate number of points
     393    Np = (m+1)*(n+1)
     394
     395    class VIndex:
     396
     397        def __init__(self, n,m):
     398            self.n = n
     399            self.m = m
     400
     401        def __call__(self, i,j):
     402            return j+i*(self.n+1)
     403
     404    class EIndex:
     405
     406        def __init__(self, n,m):
     407            self.n = n
     408            self.m = m
     409
     410        def __call__(self, i,j):
     411            return 2*(j+i*self.n)
     412
     413
     414    I = VIndex(n,m)
     415    E = EIndex(n,m)
     416
     417    points = zeros( (Np,2), Float)
     418
     419    for i in range(m+1):
     420        for j in range(n+1):
     421
     422            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
     423
     424    #Construct 2 triangles per rectangular element and assign tags to boundary
     425    #Calculate number of triangles
     426    Nt = 2*m*n
     427
     428
     429    elements = zeros( (Nt,3), Int)
     430    boundary = {}
     431    ghosts = {}
     432    nt = -1
     433    for i in range(m):
     434        for j in range(n):
     435
     436            i1 = I(i,j+1)
     437            i2 = I(i,j)
     438            i3 = I(i+1,j+1)
     439            i4 = I(i+1,j)
     440
     441            #Lower Element
     442            nt = E(i,j)
     443            if i == m-1:
     444                ghosts[nt] = E(1,j)
     445            if i == 0:
     446                ghosts[nt] = E(m-2,j)
     447
     448            if i == m-1:
     449                boundary[nt, 2] = 'right'
     450            if j == 0:
     451                boundary[nt, 1] = 'bottom'
     452            elements[nt,:] = [i4,i3,i2]
     453
     454            #Upper Element
     455            nt = E(i,j)+1
     456            if i == m-1:
     457                ghosts[nt] = E(1,j)+1
     458            if i == 0:
     459                ghosts[nt] = E(m-2,j)+1
     460
     461            if i == 0:
     462                boundary[nt, 2] = 'left'
     463            if j == n-1:
     464                boundary[nt, 1] = 'top'
     465            elements[nt,:] = [i1,i2,i3]
     466
     467
     468    return points, elements, boundary, ghosts
  • inundation/ga/storm_surge/parallel/run_parallel_advection.py

    r1407 r1414  
    1515processor_name = pypar.Get_processor_name()
    1616
     17
     18
    1719N = 30
    1820M = 30
     
    2123
    2224#Create advection domain with direction (1,-1)
    23 domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0],processor = myid)
     25domain = Parallel_Domain(points, vertices, boundary, ghosts, velocity=[1.0, 0.0])
    2426
    2527# Initial condition is zero by default
  • inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

    r1407 r1414  
    2121        self.z_models = {}
    2222        keys = self.domain.quantities.keys()
    23         print keys
     23        #print keys
    2424        for key in keys:
    2525            self.z_models[key] = faces(frame= self.frame)
    2626
    27         print self.z_models
     27        #print self.z_models
    2828
    2929        self.max_x = max(max(self.vertices[:,0]),max(self.vertices[:,2]),max(self.vertices[:,4]))
Note: See TracChangeset for help on using the changeset viewer.