Changeset 8234


Ignore:
Timestamp:
Oct 24, 2011, 5:37:21 PM (12 years ago)
Author:
steve
Message:

Working on padarn's code

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

    r8207 r8234  
    13171317        self.update_boundary()
    13181318
     1319
    13191320        # Update extrema if necessary (for reporting)
    13201321        self.update_extrema()
     
    14531454        self.update_timestep(yieldstep, finaltime)
    14541455
    1455         # Update conserved quantities
     1456        # Update centroid values of conserved quantities
    14561457        self.update_conserved_quantities()
    14571458
  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

    r8209 r8234  
    142142        #distribute_to_vertices_and_edges_limit_s_v_h(self)
    143143        distribute_to_vertices_and_edges_limit_s_v(self)
    144        
     144
     145    def update_derived_quantites(self):
     146        pass
    145147
    146148#=============== End of Channel Domain ===============================
     
    200202
    201203    w_C[:] = h_C + z_C
     204
     205    print w_C
    202206
    203207    # Extrapolate velocity and stage as well as channel geometry.
     
    217221    b_V  = width.vertex_values
    218222
    219     # These quantites need to update vertex_values
     223    # Need to update these vertex_values
    220224    a_V  = area.vertex_values
    221225    h_V  = height.vertex_values
  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

    r7884 r8234  
    66from anuga_1d.channel.channel_domain import *
    77from anuga_1d.config import g, epsilon
    8 from anuga_1d.base.generic_mesh import interval_mesh
     8from anuga_1d.base.generic_mesh import uniform_mesh
    99
    1010
     
    4242
    4343# Create domain with centroid points as defined above       
    44 domain = Domain(*interval_mesh(N))
     44domain = Domain(*uniform_mesh(N))
    4545
    4646
     
    6969
    7070
    71 exit()
     71#exit()
    7272
    7373HeightC    = domain.quantities['height'].centroid_values
  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_2_padarn.py

    r8204 r8234  
    3636import time
    3737
    38 # Set final time and yield time for simulation
    39 finaltime = 10.0
    40 yieldstep = finaltime
     38
    4139
    4240# Length of channel (m)
     
    7674t0 = time.time()
    7775
    78 finaltime= 0.0
     76# Set final time and yield time for simulation
     77finaltime = 10.0
     78yieldstep = finaltime
    7979
    8080#===================================================================
  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_depth_only.py

    r7884 r8234  
    22import random
    33from math import sqrt, pow, pi
    4 from channel_domain_Ab import *
    5 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    6 from config import g, epsilon
     4from channel_domain import *
     5from numpy import allclose, array, zeros, ones, take, sqrt
     6from anuga_1d.config import g, epsilon
    77
    88
     
    1010
    1111# Define functions for initial quantities
    12 def initial_area(x):
    13     y = zeros(len(x),Float)
    14     for i in range(len(x)):
    15         y[i]=(10-randomarray[i])
    16 
    17    
    18     return y
     12
     13
     14
     15#def initial_area(x):
     16#    y = zeros(len(x),'f')
     17#    for i in range(len(x)):
     18#        y[i]=(10-randomarray[i])
     19#    return y
     20
     21
    1922
    2023def bed(x):
    2124
    22     y = zeros(len(x),Float)
     25    y = zeros(len(x),'f')
    2326    for i in range(len(x)):
    2427        y[i]=randomarray[i]
     
    2629    return y
    2730
    28 def stage(x):
    29 
    30     y = zeros(len(x),Float)
    31     for i in range(len(x)):
    32         y[i]=3+x[i]/2000
     31
     32 
     33
     34
     35def width(x):
     36    y = zeros(len(x),'f')
     37    return y+1
     38
     39
     40stage = 6.0
     41
     42def initial_area(x):
     43
     44    a_bed = bed(x)
     45    a_width = width(x)
     46
     47    a_height = 6.0 - a_bed
     48
     49    y = a_height*a_width
     50
    3351    return y
    34  
    35 
    36 
    37 def width(x):
    38     return 1
    39  
     52
     53
    4054import time
    4155
     
    5367print "Evaluating domain with %d cells" %N
    5468cell_len = L/N # Origin = 0.0
    55 points = zeros(N+1,Float)
     69points = zeros(N+1,'f')
    5670
    5771# Define the centroid points
     
    6478# Define random array for width
    6579
    66 randomarray=zeros(len(points),Float)
     80randomarray=zeros(len(points),'f')
    6781for j in range(N+1):
    6882    randomarray[j]=random.normalvariate(3,1)
     
    7084
    7185# Set initial values of quantities - default to zero
    72 domain.set_quantity('stage',stage)
    73 #domain.set_quantity('area', initial_area)
    74 domain.set_quantity('elevation',bed)
    75 domain.set_quantity('width',width)
    76 domain.setstageflag = True
     86#domain.set_quantity('stage',6.0)
     87domain.set_quantity('elevation',bed, location='centroids')
     88domain.set_quantity('width',width, location='centroids')
     89domain.set_quantity('area', initial_area, location='centroids')
     90
     91
     92#domain.setstageflag = True
    7793# Set boundry type, order, timestepping method and limiter
    7894domain.set_boundary({'exterior':Reflective_boundary(domain)})
     
    8399#domain.h0=0.0001
    84100
     101
     102#domain.distribute_to_vertices_and_edges()
     103
     104
     105#AreaC = domain.quantities['area'].centroid_values
     106#BedC = domain.quantities['elevation'].centroid_values
     107#WidthC = domain.quantities['width'].centroid_values
     108#
     109#AreaC[:] = (6.0 - BedC)* WidthC
     110
     111#domain.set_quantity('area', initial_area)
     112
     113#domain.distribute_to_vertices_and_edges()
     114
    85115# Start timer
    86116t0 = time.time()
    87 i=1
    88 
     117i=0
     118
     119print 'elevation vertex values'
    89120print domain.quantities['elevation'].vertex_values
     121print 'stage vertex values'
     122print domain.quantities['stage'].vertex_values
     123print 'area vertex values'
     124print domain.quantities['area'].vertex_values
     125print 'width vertex values'
     126print domain.quantities['width'].vertex_values
     127
     128HeightC = domain.quantities['height'].centroid_values
     129DischargeC = domain.quantities['discharge'].centroid_values
     130C = domain.centroids
     131print 'That took %.2f seconds' %(time.time()-t0)
     132X = (domain.vertices).flatten()
     133HeightQ = (domain.quantities['height'].vertex_values).flatten()
     134VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
     135Z = (domain.quantities['elevation'].vertex_values).flatten()
     136stage=HeightQ + Z
     137
     138
     139from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
     140import pylab as p
     141#hold(False)
     142
     143plot1=subplot(211)
     144
     145plot(X,stage,X,Z)
     146
     147plot1.set_ylim([-1,11])
     148xlabel('Position')
     149ylabel('Stage')
     150legend(('Solution', 'Bed Height'),
     151                'upper right', shadow=True)
     152
     153plot2=subplot(212)
     154
     155
     156plot(X,VelocityQ)
     157plot2.set_ylim([-5,5])
     158xlabel('Position')
     159ylabel('Velocity')
     160
     161savefig(str(i)+'.png')
     162i+=1
     163plot1.clear()
     164plot2.clear()
     165
    90166
    91167for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    92168    domain.write_time()
    93169
    94     N = float(N)
    95170    HeightC = domain.quantities['height'].centroid_values
    96171    DischargeC = domain.quantities['discharge'].centroid_values
    97172    C = domain.centroids
    98173    print 'That took %.2f seconds' %(time.time()-t0)
    99     X = domain.vertices
    100     HeightQ = domain.quantities['height'].vertex_values
    101     VelocityQ = domain.quantities['velocity'].vertex_values
    102     x = X.flat
    103     z = domain.quantities['elevation'].vertex_values.flat
    104     stage=HeightQ.flat+z
    105 
    106 
     174    X = (domain.vertices).flatten()
     175    HeightQ = (domain.quantities['height'].vertex_values).flatten()
     176    VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
     177    Z = (domain.quantities['elevation'].vertex_values).flatten()
     178    stage=HeightQ + Z
    107179
    108180
     
    113185    plot1=subplot(211)
    114186
    115     plot(x,stage,x,z)
     187    plot(X,stage,X,Z)
    116188 
    117189    plot1.set_ylim([-1,11])
     
    122194
    123195    plot2=subplot(212)
    124     plot(x,VelocityQ.flat)
     196
     197   
     198    plot(X,VelocityQ)
    125199    plot2.set_ylim([-5,5])
    126200    xlabel('Position')
     
    130204    i+=1
    131205    plot1.clear()
     206    plot2.clear()
    132207
    133208print domain.quantities['elevation'].vertex_values
     
    138213C = domain.centroids
    139214print 'That took %.2f seconds' %(time.time()-t0)
    140 X = domain.vertices
    141 HeightQ = domain.quantities['height'].vertex_values
    142 VelocityQ = domain.quantities['velocity'].vertex_values
    143 x = X.flat
    144 z = domain.quantities['elevation'].vertex_values.flat
    145 stage=HeightQ.flat+z
    146 
    147 
     215X = (domain.vertices).flatten()
     216HeightQ = (domain.quantities['height'].vertex_values).flatten()
     217VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
     218Z = (domain.quantities['elevation'].vertex_values).flatten()
     219stage=HeightQ + Z
    148220
    149221
     
    154226plot1=subplot(211)
    155227
    156 plot(x,stage,x,z)
     228plot(X,stage,X,Z)
    157229 
    158230plot1.set_ylim([-1,11])
     
    163235
    164236plot2=subplot(212)
    165 plot(x,VelocityQ.flat)
     237plot(X,VelocityQ)
    166238plot2.set_ylim([-5,5])
    167239xlabel('Position')
Note: See TracChangeset for help on using the changeset viewer.