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

Continuing to numpy the for loops

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.