# Changeset 7860

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

Continuing to numpy the for loops

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

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

 r7842 self.set_timestepping_method(timestepping_method) self.wet_nodes = numpy.zeros((N,2), numpy.int) # should this be here # work array for special vertices (like wet or dry cells self.special_vertices = numpy.zeros((N,2), numpy.int) # should this be here #Allocate space for neighbour and boundary structures self.normals = numpy.zeros((N, 2), numpy.float) for i in range(N): xl = self.coordinates[i] xr = self.coordinates[i+1] self.vertices[i,0] = xl self.vertices[i,1] = xr centroid = (xl+xr)/2.0 self.centroids[i] = centroid msg = 'Coordinates should be ordered, smallest to largest' assert xr>xl, msg #The normal vectors # - point outward from each edge # - are orthogonal to the edge # - have unit length # - Are enumerated by left vertex then right vertex normals nl = -1.0 nr =  1.0 self.normals[i,:] = [nl, nr] self.areas[i] = (xr-xl) # #         print 'N', N # #         print 'Centroid', self.centroids # #         print 'Areas', self.areas # #         print 'Vertex_Coordinates', self.vertices #Initialise Neighbours (-1 means that it is a boundary neighbour) self.neighbours[i, :] = [-1, -1] #Initialise edge ids of neighbours #Initialise vertex ids of neighbours #In case of boundaries this slot is not used #self.neighbour_edges[i, :] = [-1, -1] self.neighbour_vertices[i, :] = [-1, -1] xl = self.coordinates[0:-1] xr = self.coordinates[1: ] self.vertices[:,0] = xl self.vertices[:,1] = xr self.centroids[:] = (xl+xr)/2.0 msg = 'Coordinates should be ordered, smallest to largest' assert (xr>xl).all(), msg #The normal vectors # - point outward from each edge # - are orthogonal to the edge # - have unit length # - Are enumerated by left vertex then right vertex normals self.normals[:,0] = -1.0 self.normals[:,1] =  1.0 self.areas[:] = (xr-xl) #Initialise Neighbours (-1 means that it is a boundary neighbour) self.neighbours[:, :] = -1 #Initialise edge ids of neighbours self.neighbour_vertices[:, :] = -1 #        for i in range(N): #            xl = self.coordinates[i] #            xr = self.coordinates[i+1] #            self.vertices[i,0] = xl #            self.vertices[i,1] = xr # #            centroid = (xl+xr)/2.0 #            self.centroids[i] = centroid # #            msg = 'Coordinates should be ordered, smallest to largest' #            assert xr>xl, msg # #            #The normal vectors #            # - point outward from each edge #            # - are orthogonal to the edge #            # - have unit length #            # - Are enumerated by left vertex then right vertex normals # #            nl = -1.0 #            nr =  1.0 #            self.normals[i,:] = [nl, nr] # #            self.areas[i] = (xr-xl) # ## #         print 'N', N ## #         print 'Centroid', self.centroids ## #         print 'Areas', self.areas ## #         print 'Vertex_Coordinates', self.vertices # #            #Initialise Neighbours (-1 means that it is a boundary neighbour) #            self.neighbours[i, :] = [-1, -1] #            #Initialise edge ids of neighbours #            #Initialise vertex ids of neighbours #            #In case of boundaries this slot is not used #            #self.neighbour_edges[i, :] = [-1, -1] #            self.neighbour_vertices[i, :] = [-1, -1] self.build_vertexlist() return if timestepping_method in [1, 2, 3]: self.timestepping_method = ['euler', 'rk2', 'rk3'][timestepping_method-1] return msg = '%s is an incorrect timestepping type'% timestepping_method raise Exception, msg return self.filename def get_spatial_order(self): return self.order def set_spatial_order(self, order): possible_orders = [1,2] if order in possible_orders: self.order = order return msg = '%s is an incorrect spatial order.\n'% limiter msg += 'Possible orders are: '+ ", ".join(["%s" % el for el in possible_orders]) raise Exception, msg def get_limiter(self): return self.limiter
• ## anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py

 r7855 [[-1.,  1.],[-1.,  1.],[-1.,  1.]]) def test_set_timestepping_method(self): domain = Generic_domain(self.points) domain.set_timestepping_method('euler') assert domain.timestepping_method == 'euler' domain.set_timestepping_method('rk2') assert domain.timestepping_method == 'rk2' domain.set_timestepping_method('rk3') assert domain.timestepping_method == 'rk3' domain.set_timestepping_method(1) assert domain.timestepping_method == 'euler' domain.set_timestepping_method(2) assert domain.timestepping_method == 'rk2' domain.set_timestepping_method(3) assert domain.timestepping_method == 'rk3' try: domain.set_timestepping_method(4) except: pass else: raise Exception,  'Should have raised "wrong method" error' def test_set_spatial_order(self): domain = Generic_domain(self.points) domain.set_spatial_order(1) assert domain.order == 1 domain.set_spatial_order(2) assert  domain.order == 2 try: domain.set_spatial_order(3) except: pass else: raise Exception,  'Should have raised "wrong order" error'
• ## anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

 r7855 except AssertionError: pass except: else: raise 'Should have raised "mising domain object" error'
• ## anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py

 r7855 p.sort_stats('cumulative').print_stats(30) print p
• ## anuga_work/development/2010-projects/anuga_1d/config.py

 r7839 """ epsilon = 1.0e-12 h0 = 1.0e-12 epsilon = 1.0e-6 h0 = 1.0e-6 default_boundary_tag = 'exterior' v_max = 100 #For use in domain_ext.c sound_speed = 500 max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order eta_w = 3.0e-3 #Wind stress coefficient rho_a = 1.2e-3 #Atmospheric density rho_w = 1023   #Fluid density [kg/m^3] (rho_w = 1023 for salt water) #Betas [0;1] control the allowed steepness of gradient for second order #extrapolations. Values of 1 allow the steepes gradients while #lower values are more conservative. Values of 0 correspond to #1'st order extrapolations. # # Large values of beta_h may cause simulations to require more timesteps # as surface will 'hug' closer to the bed. # Small values of beta_h will make code faster, but one may experience # artificial momenta caused by discontinuities in water depths in # the presence of steep slopes. One example of this would be # stationary water 'lapping' upwards to a higher point on the coast. # # # #There are separate betas for the w-limiter and the h-limiter # # # # #Good values are: #beta_w = 0.9 #beta_h = 0.2 beta_w = 1.5 beta_h = 0.2 timestepping_method = 'euler' beta_w = 1.0 beta_h = 1.0 timestepping_method = 'rk2' CFL = 1.0  #FIXME (ole): Is this in use yet?? use_extensions = True    #Try to use C-extensions #use_extensions = False   #Do not use C-extensions use_psyco = True  #Use psyco optimisations #use_psyco = False  #Do not use psyco optimisations optimised_gradient_limiter = True #Use hardwired gradient limiter #use_extensions = True    #Try to use C-extensions ##use_extensions = False   #Do not use C-extensions # #use_psyco = True  #Use psyco optimisations ##use_psyco = False  #Do not use psyco optimisations # # #optimised_gradient_limiter = True #Use hardwired gradient limiter #Specific to shallow water W.E.
• ## anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

 r7857 import os from math import sqrt, pi from shallow_water_vel_domain import * from Numeric import allclose, array, zeros, ones, Float, take, sqrt from config import g, epsilon import numpy import time #from Numeric import allclose, array, zeros, ones, Float, take, sqrt from anuga_1d.sww.sww_domain import * from anuga_1d.config import g, epsilon from anuga_1d.base.generic_mesh import uniform_mesh h1 = 10.0 return h , u*h, u print "TEST 1D-SOLUTION III -- DRY BED" def stage(x): y = zeros(len(x),Float) for i in range(len(x)): if x[i]<=L/4.0: y[i] = h0 elif x[i]<=3*L/4.0: y[i] = h1 else: y[i] = h0 import numpy y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0) #    for i in range(len(x)): #        if x[i]<=L/4.0: #            y[i] = h0 #        elif x[i]<=3*L/4.0: #            y[i] = h1 #        else: #            y[i] = h0 return y import time finaltime = 2.0 print "TEST 1D-SOLUTION III -- DRY BED" finaltime = 4.0 yieldstep = 0.1 L = 2000.0     # Length of channel (m) N = 800 print "Evaluating domain with %d cells" %N cell_len = L/N # Origin = 0.0 points = zeros(N+1,Float) for j in range(N+1): points[j] = j*cell_len domain = Domain(points) domain = Domain(*uniform_mesh(N)) domain.set_quantity('stage', stage) domain.set_boundary({'exterior': Reflective_boundary(domain)}) Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right' : Br}) domain.order = 2 domain.set_timestepping_method('euler') domain.set_CFL(1.0) domain.set_limiter("vanleer") domain.set_limiter("minmod") #domain.h0=0.0001 domain.write_time() print 'end' print 'That took %.2f seconds' %(time.time()-t0) N = float(N) HeightC    = domain.quantities['height'].centroid_values StageC     = domain.quantities['stage'].centroid_values BedC       = domain.quantities['elevation'].centroid_values C          = domain.centroids HeightV    = domain.quantities['height'].vertex_values StageV     = domain.quantities['stage'].vertex_values BedV       = domain.quantities['elevation'].vertex_values VelocityV  = domain.quantities['velocity'].vertex_values X          = domain.vertices import pylab pylab.hold(False) plot1 = pylab.subplot(211) pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat) plot1.set_ylim([-1,11]) pylab.xlabel('Position') pylab.ylabel('Stage') pylab.legend(('Analytical Solution', 'Numerical Solution'), 'upper right', shadow=True) plot2 = pylab.subplot(212) pylab.plot(X.flat,VelocityV.flat) plot2.set_ylim([-20,10]) pylab.xlabel('Position') pylab.ylabel('Velocity') pylab.show()
• ## anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py

 r7852 def gravity(domain): def gravity_for_loops(domain): """Apply gravitational pull in the presence of bed slope """ xmom  = domain.quantities['xmomentum'].explicit_update stage = domain.quantities['stage'].explicit_update #    ymom = domain.quantities['ymomentum'].explicit_update Stage = domain.quantities['stage'] Elevation = domain.quantities['elevation'] #h = Stage.edge_values - Elevation.edge_values h = Stage.vertex_values - Elevation.vertex_values b = Elevation.vertex_values for k in range(domain.number_of_elements): #        avg_h = sum( h[k,:] )/3 avg_h = sum( h[k,:] )/2 b0, b1 = b[k,:] w0, w1 = w[k,:] wx = gradient(x0, x1, w0, w1) #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) #stage[k] = 0.0 def gravity(domain): """Apply gravitational pull in the presence of bed slope """ xmom  = domain.quantities['xmomentum'].explicit_update Elevation = domain.quantities['elevation'] Height    = domain.quantities['height'] hc = Height.centroid_values b  = Elevation.vertex_values x = domain.vertices g = domain.g x0 = x[:,0] x1 = x[:,1] b0 = b[:,0] b1 = b[:,1] bx = (b1-b0)/(x1-x0) xmom[:] = xmom - g*bx*hc def manning_friction(domain):
• ## anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_domain.py

 r7852 w_C[i] = z_C[i] else: #u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i]) u_C[i]  = uh_C[i]/h_C[i] u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i]) #u_C[i]  = uh_C[i]/h_C[i] for name in [ 'velocity', 'stage' ]:
• ## anuga_work/development/2010-projects/anuga_1d/sww/test_sww_domain.py

 r7839 def test_gravity(self): def test_gravity_for_loops(self): """ Compare shallow_water_domain gravity calculation domain.set_boundary({'exterior' : Reflective_boundary(domain)}) domain.set_spatial_order(1) domain.distribute_to_vertices_and_edges() gravity_for_loops(domain) assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update ) def test_gravity(self): """ Compare shallow_water_domain gravity calculation """ def slope_one(x): return x domain = Domain(self.points) domain.set_quantity('stage',4.0) domain.set_quantity('elevation',slope_one) domain.set_boundary({'exterior' : Reflective_boundary(domain)}) domain.set_spatial_order(1) domain.distribute_to_vertices_and_edges() gravity(domain) #print domain.quantities['stage'].vertex_values #print domain.quantities['elevation'].vertex_values #print domain.quantities['xmomentum'].explicit_update assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update ) if __name__ == "__main__": suite = unittest.makeSuite(Test_Shallow_Water, 'test') suite = unittest.makeSuite(Test_Shallow_Water, 'test_gravity') #suite = unittest.makeSuite(Test_Shallow_Water, 'test_evolve_first_order')
• ## anuga_work/development/2010-projects/anuga_1d/test_all.py

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