Changeset 7860


Ignore:
Timestamp:
Jun 18, 2010, 5:49:57 PM (15 years ago)
Author:
steve
Message:

Continuing to numpy the for loops

Location:
anuga_work/development/2010-projects/anuga_1d
Files:
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py

    r7842 r7860  
    4040        self.set_timestepping_method(timestepping_method)
    4141
    42         self.wet_nodes = numpy.zeros((N,2), numpy.int) # should this be here
     42
     43        # work array for special vertices (like wet or dry cells
     44        self.special_vertices = numpy.zeros((N,2), numpy.int) # should this be here
    4345
    4446        #Allocate space for neighbour and boundary structures
    … …  
    5860
    5961        self.normals = numpy.zeros((N, 2), numpy.float)
    60        
    61         for i in range(N):
    62             xl = self.coordinates[i]
    63             xr = self.coordinates[i+1]
    64             self.vertices[i,0] = xl
    65             self.vertices[i,1] = xr
    66 
    67             centroid = (xl+xr)/2.0
    68             self.centroids[i] = centroid
    69 
    70             msg = 'Coordinates should be ordered, smallest to largest'
    71             assert xr>xl, msg
    72            
    73             #The normal vectors
    74             # - point outward from each edge
    75             # - are orthogonal to the edge
    76             # - have unit length
    77             # - Are enumerated by left vertex then right vertex normals
    78 
    79             nl = -1.0
    80             nr =  1.0
    81             self.normals[i,:] = [nl, nr]
    82 
    83             self.areas[i] = (xr-xl)
    84            
    85 # #         print 'N', N
    86 # #         print 'Centroid', self.centroids
    87 # #         print 'Areas', self.areas
    88 # #         print 'Vertex_Coordinates', self.vertices
    89            
    90             #Initialise Neighbours (-1 means that it is a boundary neighbour)
    91             self.neighbours[i, :] = [-1, -1]
    92             #Initialise edge ids of neighbours
    93             #Initialise vertex ids of neighbours
    94             #In case of boundaries this slot is not used
    95             #self.neighbour_edges[i, :] = [-1, -1]
    96             self.neighbour_vertices[i, :] = [-1, -1]
     62
     63
     64        xl = self.coordinates[0:-1]
     65        xr = self.coordinates[1: ]
     66
     67        self.vertices[:,0] = xl
     68        self.vertices[:,1] = xr
     69
     70        self.centroids[:] = (xl+xr)/2.0
     71
     72        msg = 'Coordinates should be ordered, smallest to largest'
     73        assert (xr>xl).all(), msg
     74
     75        #The normal vectors
     76        # - point outward from each edge
     77        # - are orthogonal to the edge
     78        # - have unit length
     79        # - Are enumerated by left vertex then right vertex normals
     80
     81
     82        self.normals[:,0] = -1.0
     83        self.normals[:,1] =  1.0
     84
     85        self.areas[:] = (xr-xl)
     86
     87        #Initialise Neighbours (-1 means that it is a boundary neighbour)
     88
     89        self.neighbours[:, :] = -1
     90
     91        #Initialise edge ids of neighbours
     92        self.neighbour_vertices[:, :] = -1
     93
     94
     95#        for i in range(N):
     96#            xl = self.coordinates[i]
     97#            xr = self.coordinates[i+1]
     98#            self.vertices[i,0] = xl
     99#            self.vertices[i,1] = xr
     100#
     101#            centroid = (xl+xr)/2.0
     102#            self.centroids[i] = centroid
     103#
     104#            msg = 'Coordinates should be ordered, smallest to largest'
     105#            assert xr>xl, msg
     106#
     107#            #The normal vectors
     108#            # - point outward from each edge
     109#            # - are orthogonal to the edge
     110#            # - have unit length
     111#            # - Are enumerated by left vertex then right vertex normals
     112#
     113#            nl = -1.0
     114#            nr =  1.0
     115#            self.normals[i,:] = [nl, nr]
     116#
     117#            self.areas[i] = (xr-xl)
     118#
     119## #         print 'N', N
     120## #         print 'Centroid', self.centroids
     121## #         print 'Areas', self.areas
     122## #         print 'Vertex_Coordinates', self.vertices
     123#
     124#            #Initialise Neighbours (-1 means that it is a boundary neighbour)
     125#            self.neighbours[i, :] = [-1, -1]
     126#            #Initialise edge ids of neighbours
     127#            #Initialise vertex ids of neighbours
     128#            #In case of boundaries this slot is not used
     129#            #self.neighbour_edges[i, :] = [-1, -1]
     130#            self.neighbour_vertices[i, :] = [-1, -1]
    97131
    98132        self.build_vertexlist()
    … …  
    630664            return
    631665
     666        if timestepping_method in [1, 2, 3]:
     667            self.timestepping_method = ['euler', 'rk2', 'rk3'][timestepping_method-1]
     668            return
     669
    632670        msg = '%s is an incorrect timestepping type'% timestepping_method
    633671        raise Exception, msg
    … …  
    814852        return self.filename
    815853
     854
     855    def get_spatial_order(self):
     856        return self.order
     857
     858    def set_spatial_order(self, order):
     859
     860        possible_orders = [1,2]
     861
     862        if order in possible_orders:
     863            self.order = order
     864            return
     865
     866        msg = '%s is an incorrect spatial order.\n'% limiter
     867        msg += 'Possible orders are: '+ ", ".join(["%s" % el for el in possible_orders])
     868        raise Exception, msg
     869   
     870
    816871    def get_limiter(self):
    817872        return self.limiter

  • anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py

    r7855 r7860  
    4646            [[-1.,  1.],[-1.,  1.],[-1.,  1.]])
    4747
     48    def test_set_timestepping_method(self):
     49
     50        domain = Generic_domain(self.points)
     51
     52        domain.set_timestepping_method('euler')
     53        assert domain.timestepping_method == 'euler'
     54
     55        domain.set_timestepping_method('rk2')
     56        assert domain.timestepping_method == 'rk2'
     57
     58        domain.set_timestepping_method('rk3')
     59        assert domain.timestepping_method == 'rk3'
     60
     61        domain.set_timestepping_method(1)
     62        assert domain.timestepping_method == 'euler'
     63
     64        domain.set_timestepping_method(2)
     65        assert domain.timestepping_method == 'rk2'
     66
     67        domain.set_timestepping_method(3)
     68        assert domain.timestepping_method == 'rk3'
     69
     70        try:
     71            domain.set_timestepping_method(4)
     72        except:
     73            pass
     74        else:
     75            raise Exception,  'Should have raised "wrong method" error'
     76
     77    def test_set_spatial_order(self):
     78
     79        domain = Generic_domain(self.points)
     80
     81        domain.set_spatial_order(1)
     82        assert domain.order == 1
     83
     84        domain.set_spatial_order(2)
     85        assert  domain.order == 2
     86
     87
     88        try:
     89            domain.set_spatial_order(3)
     90        except:
     91            pass
     92        else:
     93            raise Exception,  'Should have raised "wrong order" error'
    4894
    4995

  • anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

    r7855 r7860  
    6565        except AssertionError:
    6666            pass
    67         except:
     67        else:
    6868            raise 'Should have raised "mising domain object" error'
    6969

  • anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py

    r7855 r7860  
    9494p.sort_stats('cumulative').print_stats(30)
    9595
    96 print p
    9796
    9897

  • anuga_work/development/2010-projects/anuga_1d/config.py

    r7839 r7860  
    22"""
    33
    4 epsilon = 1.0e-12
    5 h0 = 1.0e-12
     4epsilon = 1.0e-6
     5h0 = 1.0e-6
    66
    77default_boundary_tag = 'exterior'
    … …  
    3333
    3434
    35 v_max = 100 #For use in domain_ext.c
    36 sound_speed = 500
    37 
    3835
    3936max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order
    … …  
    4744
    4845
    49 eta_w = 3.0e-3 #Wind stress coefficient
    50 rho_a = 1.2e-3 #Atmospheric density
    51 rho_w = 1023   #Fluid density [kg/m^3] (rho_w = 1023 for salt water)
    5246
    5347
    54 #Betas [0;1] control the allowed steepness of gradient for second order
    55 #extrapolations. Values of 1 allow the steepes gradients while
    56 #lower values are more conservative. Values of 0 correspond to
    57 #1'st order extrapolations.
    58 #
    59 # Large values of beta_h may cause simulations to require more timesteps
    60 # as surface will 'hug' closer to the bed.
    61 # Small values of beta_h will make code faster, but one may experience
    62 # artificial momenta caused by discontinuities in water depths in
    63 # the presence of steep slopes. One example of this would be
    64 # stationary water 'lapping' upwards to a higher point on the coast.
    65 #
    66 #
    67 #
    68 #There are separate betas for the w-limiter and the h-limiter
    69 #
    70 #
    71 #
    72 #
    73 #Good values are:
    74 #beta_w = 0.9
    75 #beta_h = 0.2
    76 
    77 
    78 
    79 beta_w = 1.5
    80 beta_h = 0.2
    81 timestepping_method = 'euler'
     48beta_w = 1.0
     49beta_h = 1.0
     50timestepping_method = 'rk2'
    8251
    8352CFL = 1.0  #FIXME (ole): Is this in use yet??
    … …  
    9968
    10069
    101 use_extensions = True    #Try to use C-extensions
    102 #use_extensions = False   #Do not use C-extensions
    103 
    104 use_psyco = True  #Use psyco optimisations
    105 #use_psyco = False  #Do not use psyco optimisations
    106 
    107 
    108 optimised_gradient_limiter = True #Use hardwired gradient limiter
     70#use_extensions = True    #Try to use C-extensions
     71##use_extensions = False   #Do not use C-extensions
     72#
     73#use_psyco = True  #Use psyco optimisations
     74##use_psyco = False  #Do not use psyco optimisations
     75#
     76#
     77#optimised_gradient_limiter = True #Use hardwired gradient limiter
    10978
    11079#Specific to shallow water W.E.

  • anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

    r7857 r7860  
    11import os
    22from math import sqrt, pi
    3 from shallow_water_vel_domain import *
    4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    5 from config import g, epsilon
     3import numpy
     4import time
     5#from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    66
     7
     8
     9from anuga_1d.sww.sww_domain import *
     10from anuga_1d.config import g, epsilon
     11from anuga_1d.base.generic_mesh import uniform_mesh
    712
    813h1 = 10.0
    … …  
    4954    return h , u*h, u
    5055
    51 print "TEST 1D-SOLUTION III -- DRY BED"
     56
    5257
    5358def stage(x):
    54     y = zeros(len(x),Float)
    55     for i in range(len(x)):
    56         if x[i]<=L/4.0:
    57             y[i] = h0
    58         elif x[i]<=3*L/4.0:
    59             y[i] = h1
    60         else:
    61             y[i] = h0
     59    import numpy
     60
     61    y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0)
     62
     63#    for i in range(len(x)):
     64#        if x[i]<=L/4.0:
     65#            y[i] = h0
     66#        elif x[i]<=3*L/4.0:
     67#            y[i] = h1
     68#        else:
     69#            y[i] = h0
    6270    return y
    6371
    6472
    65 import time
    6673
    67 finaltime = 2.0
     74print "TEST 1D-SOLUTION III -- DRY BED"
     75
     76
     77finaltime = 4.0
    6878yieldstep = 0.1
    6979L = 2000.0     # Length of channel (m)
    … …  
    7383N = 800
    7484print "Evaluating domain with %d cells" %N
    75 cell_len = L/N # Origin = 0.0
    76 points = zeros(N+1,Float)
    77 
    78 for j in range(N+1):
    79     points[j] = j*cell_len
    80        
    81 domain = Domain(points)
     85domain = Domain(*uniform_mesh(N))
    8286   
    8387domain.set_quantity('stage', stage)
    84 domain.set_boundary({'exterior': Reflective_boundary(domain)})
     88
     89Br = Reflective_boundary(domain)
     90
     91domain.set_boundary({'left': Br, 'right' : Br})
    8592domain.order = 2
    8693domain.set_timestepping_method('euler')
    8794domain.set_CFL(1.0)
    88 domain.set_limiter("vanleer")
     95domain.set_limiter("minmod")
    8996#domain.h0=0.0001
    9097
    … …  
    94101    domain.write_time()
    95102
    96 print 'end'
     103print 'That took %.2f seconds' %(time.time()-t0)
    97104
    98105
     106N = float(N)
     107HeightC    = domain.quantities['height'].centroid_values
     108StageC     = domain.quantities['stage'].centroid_values
     109BedC       = domain.quantities['elevation'].centroid_values
     110C          = domain.centroids
     111
     112HeightV    = domain.quantities['height'].vertex_values
     113StageV     = domain.quantities['stage'].vertex_values
     114BedV       = domain.quantities['elevation'].vertex_values
     115VelocityV  = domain.quantities['velocity'].vertex_values
     116X          = domain.vertices
     117
     118
     119import pylab
     120
     121pylab.hold(False)
     122
     123plot1 = pylab.subplot(211)
     124
     125pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
     126
     127plot1.set_ylim([-1,11])
     128
     129pylab.xlabel('Position')
     130pylab.ylabel('Stage')
     131pylab.legend(('Analytical Solution', 'Numerical Solution'),
     132           'upper right', shadow=True)
     133
     134
     135plot2 = pylab.subplot(212)
     136
     137pylab.plot(X.flat,VelocityV.flat)
     138plot2.set_ylim([-20,10])
     139
     140pylab.xlabel('Position')
     141pylab.ylabel('Velocity')
     142
     143pylab.show()
     144

  • anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py

    r7852 r7860  
    1010
    1111
    12 def gravity(domain):
     12def gravity_for_loops(domain):
    1313    """Apply gravitational pull in the presence of bed slope
    1414    """
    … …  
    1717
    1818    xmom  = domain.quantities['xmomentum'].explicit_update
    19     stage = domain.quantities['stage'].explicit_update
    20 #    ymom = domain.quantities['ymomentum'].explicit_update
    2119
    2220    Stage = domain.quantities['stage']
    2321    Elevation = domain.quantities['elevation']
    24     #h = Stage.edge_values - Elevation.edge_values
     22
    2523    h = Stage.vertex_values - Elevation.vertex_values
    2624    b = Elevation.vertex_values
    … …  
    3129
    3230    for k in range(domain.number_of_elements):
    33 #        avg_h = sum( h[k,:] )/3
    3431        avg_h = sum( h[k,:] )/2
    3532
    … …  
    4037        b0, b1 = b[k,:]
    4138
    42         w0, w1 = w[k,:]
    43         wx = gradient(x0, x1, w0, w1)
    4439
    4540        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
    … …  
    5146        #stage[k] = 0.0
    5247
     48
     49def gravity(domain):
     50    """Apply gravitational pull in the presence of bed slope
     51    """
     52
     53    xmom  = domain.quantities['xmomentum'].explicit_update
     54
     55    Elevation = domain.quantities['elevation']
     56    Height    = domain.quantities['height']
     57
     58    hc = Height.centroid_values
     59    b  = Elevation.vertex_values
     60
     61    x = domain.vertices
     62    g = domain.g
     63
     64    x0 = x[:,0]
     65    x1 = x[:,1]
     66
     67    b0 = b[:,0]
     68    b1 = b[:,1]
     69
     70    bx = (b1-b0)/(x1-x0)
     71
     72    xmom[:] = xmom - g*bx*hc
    5373
    5474def manning_friction(domain):

  • anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_domain.py

    r7852 r7860  
    246246            w_C[i] = z_C[i]
    247247        else:
    248             #u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
    249             u_C[i]  = uh_C[i]/h_C[i]
     248            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
     249            #u_C[i]  = uh_C[i]/h_C[i]
    250250       
    251251    for name in [ 'velocity', 'stage' ]:

  • anuga_work/development/2010-projects/anuga_1d/sww/test_sww_domain.py

    r7839 r7860  
    7777
    7878
    79     def test_gravity(self):
     79    def test_gravity_for_loops(self):
    8080        """
    8181        Compare shallow_water_domain gravity calculation
    … …  
    9090        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    9191
     92        domain.set_spatial_order(1)
     93        domain.distribute_to_vertices_and_edges()
     94
     95        gravity_for_loops(domain)     
     96
     97        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
     98
     99
     100    def test_gravity(self):
     101        """
     102        Compare shallow_water_domain gravity calculation
     103        """
     104
     105        def slope_one(x):
     106            return x
     107
     108        domain = Domain(self.points)
     109        domain.set_quantity('stage',4.0)
     110        domain.set_quantity('elevation',slope_one)
     111        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
     112
     113
     114        domain.set_spatial_order(1)
     115        domain.distribute_to_vertices_and_edges()
     116
    92117        gravity(domain)
    93        
    94         #print domain.quantities['stage'].vertex_values
    95         #print domain.quantities['elevation'].vertex_values
    96         #print domain.quantities['xmomentum'].explicit_update       
    97118
    98119        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )
    … …  
    180201
    181202if __name__ == "__main__":
    182     suite = unittest.makeSuite(Test_Shallow_Water, 'test')
     203    suite = unittest.makeSuite(Test_Shallow_Water, 'test_gravity')
    183204    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
    184205

  • anuga_work/development/2010-projects/anuga_1d/test_all.py

    r7855 r7860  
    191191    check_anuga_1d_import()
    192192
    193     print 'test-all'
    194 
    195193    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
    196194        test_verbose = True
Note: See TracChangeset for help on using the changeset viewer.