Changeset 8241


Ignore:
Timestamp:
Oct 27, 2011, 8:23:34 PM (14 years ago)
Author:
steve
Message:

Got the ctac figures working

Location:
trunk/anuga_work/development/2010-projects/anuga_1d/channel
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

    r8240 r8241  
    151151    def initialize_plotting(self,
    152152                            stage_lim = [-1.0, 40.0],
    153                             height_lim = [-1.0, 10.0],
     153                            width_lim = [-1.0, 10.0],
    154154                            velocity_lim = [-5.0, 5.0]):
    155155
    … …  
    159159
    160160        x = self.get_vertices().flatten()
     161
    161162        z = self.quantities['elevation'].vertex_values.flatten()
    162163        w = self.quantities['stage'].vertex_values.flatten()
    163164        h = self.quantities['height'].vertex_values.flatten()
    164165        v = self.quantities['velocity'].vertex_values.flatten()
     166        b = self.quantities['width'].vertex_values.flatten()
    165167
    166168        print x.shape
    167169        print z.shape
    168170
     171        #-------------------------------
     172        # Top plot
     173        #-------------------------------
    169174        self.plot1 = pylab.subplot(311)
    170175
    … …  
    173178
    174179        self.plot1.set_ylim(stage_lim)
    175         pylab.xlabel('Position')
     180        #pylab.xlabel('Position')
    176181        pylab.ylabel('Stage')
    177182
     183        #-------------------------------
     184        # Middle Plot
     185        #-------------------------------
    178186        self.plot2 = pylab.subplot(312)
    179187
    180         self.hplot, = pylab.plot(x, h)
    181 
    182         self.plot2.set_ylim(height_lim)
    183         pylab.xlabel('Position')
    184         pylab.ylabel('Height')
    185 
     188        self.bplot, = pylab.plot(x, b)
     189
     190        self.plot2.set_ylim(width_lim)
     191        #pylab.xlabel('Position')
     192        pylab.ylabel('Width')
     193
     194        #-------------------------------
     195        # Bottom Plot
     196        #-------------------------------
    186197        self.plot3 = pylab.subplot(313)
    187198
    … …  
    203214        h = self.quantities['height'].vertex_values.flatten()
    204215        v = self.quantities['velocity'].vertex_values.flatten()
     216        b = self.quantities['width'].vertex_values.flatten()
    205217
    206218
    207219        self.zplot.set_ydata(z)
    208220        self.wplot.set_ydata(w)
    209         self.hplot.set_ydata(h)
     221        self.bplot.set_ydata(b)
    210222        self.vplot.set_ydata(v)
    211223
    … …  
    213225
    214226
    215     def hold_plotting(self):
     227    def hold_plotting(self,save=None):
    216228
    217229        self.update_plotting()
    … …  
    219231       
    220232        pylab.ioff()
     233
     234        if save != None:
     235            file = save+".pdf"
     236            pylab.savefig(file)
    221237
    222238        pylab.show()
    … …  
    287303    w_C[:] = h_C + z_C
    288304
    289     print w_C
     305    #print w_C
    290306
    291307    # Extrapolate velocity and stage as well as channel geometry.

  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_domain_ext.c

    r7884 r8241  
    624624            qrm[5] = width_edge_values[ki];
    625625
    626              _flux_function_channel(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
     626             _flux_function_channel21(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
    627627            flux[0] -= edgeflux[0];
    628628            flux[1] -= edgeflux[1];

  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_depth_only.py

    r8240 r8241  
    105105WidthC = domain.quantities['width'].centroid_values
    106106#
    107 AreaC[:] = (6.0 - BedC)* WidthC
     107AreaC[:] = (8.0 - BedC)* WidthC
    108108
    109109#domain.set_quantity('area', initial_area)
    … …  
    124124print domain.quantities['width'].vertex_values
    125125
    126 #HeightC = domain.quantities['height'].centroid_values
    127 #DischargeC = domain.quantities['discharge'].centroid_values
    128 #C = domain.centroids
    129 #print 'That took %.2f seconds' %(time.time()-t0)
    130 #X = (domain.vertices).flatten()
    131 #HeightQ = (domain.quantities['height'].vertex_values).flatten()
    132 #VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
    133 #Z = (domain.quantities['elevation'].vertex_values).flatten()
    134 #stage=HeightQ + Z
    135 #
    136 #
    137 #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
    138 #import pylab as p
    139 ##hold(False)
    140 #
    141 #plot1=subplot(211)
    142 #
    143 #plot(X,stage,X,Z)
    144 #
    145 #plot1.set_ylim([-1,11])
    146 #xlabel('Position')
    147 #ylabel('Stage')
    148 #legend(('Solution', 'Bed Height'),
    149 #                'upper right', shadow=True)
    150 #
    151 #plot2=subplot(212)
    152 #
    153 #
    154 #plot(X,VelocityQ)
    155 #plot2.set_ylim([-5,5])
    156 #xlabel('Position')
    157 #ylabel('Velocity')
    158 #
    159 #savefig(str(i)+'.png')
    160 #i+=1
    161 #plot1.clear()
    162 #plot2.clear()
    163126
    164127domain.distribute_to_vertices_and_edges()
    … …  
    170133yieldstep = 1.0
    171134
    172 
    173 
    174 domain.initialize_plotting(stage_lim = [-2.0, 14.0],
    175                            height_lim = [-2.0, 14.0],
     135domain.initialize_plotting(stage_lim = [-2.0, 12.0],
     136                           width_lim = [0.0, 2.0],
    176137                           velocity_lim = [-10.0, 10.0])
    177138
    … …  
    181142    domain.update_plotting()
    182143
    183 #    HeightC = domain.quantities['height'].centroid_values
    184 #    DischargeC = domain.quantities['discharge'].centroid_values
    185 #    C = domain.centroids
    186 #    print 'That took %.2f seconds' %(time.time()-t0)
    187 #    X = (domain.vertices).flatten()
    188 #    HeightQ = (domain.quantities['height'].vertex_values).flatten()
    189 #    VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
    190 #    Z = (domain.quantities['elevation'].vertex_values).flatten()
    191 #    stage=HeightQ + Z
    192 #
    193 #
    194 #    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
    195 #    import pylab as p
    196 ##hold(False)
    197 #
    198 #    plot1=subplot(211)
    199 #
    200 #    plot(X,stage,X,Z)
    201 #
    202 #    plot1.set_ylim([-1,11])
    203 #    xlabel('Position')
    204 #    ylabel('Stage')
    205 #    legend(('Solution', 'Bed Height'),
    206 #                'upper right', shadow=True)
    207 #
    208 #    plot2=subplot(212)
    209 #
    210 #
    211 #    plot(X,VelocityQ)
    212 #    plot2.set_ylim([-5,5])
    213 #    xlabel('Position')
    214 #    ylabel('Velocity')
    215 #
    216 #    savefig(str(i)+'.png')
    217 #    i+=1
    218 #    plot1.clear()
    219 #    plot2.clear()
    220144
    221145print domain.quantities['elevation'].vertex_values
    222146
    223 
    224 domain.hold_plotting()
     147domain.hold_plotting(save="not_well_balanced_random_depth")
    225148
    226149
    227 #N = float(N)
    228 #HeightC = domain.quantities['height'].centroid_values
    229 #DischargeC = domain.quantities['discharge'].centroid_values
    230 #C = domain.centroids
    231 #print 'That took %.2f seconds' %(time.time()-t0)
    232 #X = (domain.vertices).flatten()
    233 #HeightQ = (domain.quantities['height'].vertex_values).flatten()
    234 #VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
    235 #Z = (domain.quantities['elevation'].vertex_values).flatten()
    236 #stage=HeightQ + Z
    237 #
    238 #
    239 #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    240 #
    241 ##hold(False)
    242 #
    243 #plot1=subplot(211)
    244 #
    245 #plot(X,stage,X,Z)
    246 #
    247 #plot1.set_ylim([-1,11])
    248 #xlabel('Position')
    249 #ylabel('Stage')
    250 #legend(('Solution', 'Bed Height'),
    251 #           'upper right', shadow=True)
    252 #
    253 #plot2=subplot(212)
    254 #plot(X,VelocityQ)
    255 #plot2.set_ylim([-5,5])
    256 #xlabel('Position')
    257 #ylabel('Velocity')
    258    
    259 #show()
    260 

  • trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_width_depth.py

    r7884 r8241  
    11import os
    22import random
    3 from 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
     3from math import sqrt, pow, pi
     4from channel_domain import *
     5from numpy import allclose, array, zeros, ones, take, sqrt
     6from anuga_1d.config import g, epsilon
    77
     8
     9random.seed(111)
    810
    911print "Variable Width Only Test"
    … …  
    1113# Define functions for initial quantities
    1214def initial_area(x):
    13     y = zeros(len(x),Float)
     15    y = zeros(len(x),'f')
    1416    for i in range(len(x)):
    1517        y[i]=(10-randomarray[i])
    … …  
    2022def width(x):
    2123
    22     y = zeros(len(x),Float)
     24    y = zeros(len(x),'f')
    2325    for i in range(len(x)):
    2426        y[i]=randomarray[i]
    … …  
    2830def bed(x):
    2931
    30     y = zeros(len(x),Float)
     32    y = zeros(len(x),'f')
    3133    for i in range(len(x)):
    3234        y[i]=randomarray2[i]
    … …  
    3840def stage(x):
    3941
    40     y = zeros(len(x),Float)
     42    y = zeros(len(x),'f')
    4143    for i in range(len(x)):
    4244        if x[i]<200 and x[i]>150:
    … …  
    6466print "Evaluating domain with %d cells" %N
    6567cell_len = L/N # Origin = 0.0
    66 points = zeros(N+1,Float)
     68points = zeros(N+1,'f')
    6769
    6870# Define the centroid points
    … …  
    7476
    7577# Define random array for width
    76 randomarray=zeros(len(points),Float)
     78randomarray=zeros(len(points),'f')
    7779for j in range(N+1):
    7880    randomarray[j]=random.normalvariate(3,1)
    79 randomarray2=zeros(len(points),Float)
     81randomarray2=zeros(len(points),'f')
    8082for j in range(N+1):
    8183    randomarray2[j]=random.normalvariate(3,1)
    … …  
    98100t0 = time.time()
    99101
    100 print domain.quantities['elevation'].vertex_values
     102
     103AreaC = domain.quantities['area'].centroid_values
     104BedC = domain.quantities['elevation'].centroid_values
     105WidthC = domain.quantities['width'].centroid_values
     106#
     107AreaC[:] = (10.0 - BedC)* WidthC
     108
     109
     110#print domain.quantities['elevation'].vertex_values
     111
     112finaltime = 100.0
     113yieldstep = 1.0
     114
     115domain.initialize_plotting(stage_lim = [-1.0, 20.0],
     116                           width_lim = [0.0, 6.0],
     117                           velocity_lim = [-1.0, 1.0])
     118
    101119
    102120for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    103121    domain.write_time()
    104122
    105 print domain.quantities['elevation'].vertex_values
     123    domain.update_plotting()
    106124
    107 N = float(N)
    108 HeightC = domain.quantities['height'].centroid_values
    109 DischargeC = domain.quantities['discharge'].centroid_values
    110 C = domain.centroids
    111 print 'That took %.2f seconds' %(time.time()-t0)
    112 X = domain.vertices
    113 HeightQ = domain.quantities['height'].vertex_values
    114 VelocityQ = domain.quantities['velocity'].vertex_values
    115 x = X.flat
    116 z = domain.quantities['elevation'].vertex_values.flat
    117 b = domain.quantities['width'].vertex_values.flat
    118 stage=HeightQ.flat+z
     125#print domain.quantities['elevation'].vertex_values
    119126
     127domain.hold_plotting()
    120128
    121 
    122 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    123 
    124 #hold(False)
    125    
    126 plot1=subplot(211)
    127 
    128 plot(x,stage,x,z,x,b)
    129  
    130 plot1.set_ylim([-5,15])
    131 xlabel('Position')
    132 ylabel('Stage')
    133 legend(('Solution', 'Bed','Width'),
    134             'upper right', shadow=True)
    135 
    136 plot2=subplot(212)
    137 plot(x,VelocityQ.flat)
    138 plot2.set_ylim([-10,10])
    139 xlabel('Position')
    140 ylabel('Velocity')
    141    
    142 show()
    143 
Note: See TracChangeset for help on using the changeset viewer.