Changeset 8265


Ignore:
Timestamp:
Dec 6, 2011, 11:09:34 AM (12 years ago)
Author:
paul
Message:

Added mass conservation check to sqpipe

Location:
trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py

    r8263 r8265  
    100100    v = domain.get_quantity('velocity', 'centroids')
    101101    t = domain.get_quantity('top', 'centroids')
    102 
    103     plot1 = pylab.subplot(311)
    104 
     102    s = domain.state
     103    m = domain.get_mass()
     104    M = m.sum() * numpy.ones_like(x)
     105
     106    plot1 = pylab.subplot(511)
    105107    zplot, = pylab.plot(x, z)
    106108    wplot, = pylab.plot(x, w)
    107     ztplot, = pylab.plot(x, z+t)
    108 
     109    ztplot, = pylab.plot(x, z+t)   
    109110    plot1.set_ylim([-1,50])
    110111    pylab.xlabel('Position')
    111112    pylab.ylabel('Stage')
    112113
    113     plot2 = pylab.subplot(312)
    114 
     114    plot2 = pylab.subplot(512)
    115115    hplot, = pylab.plot(x, h)
    116116    tplot, = pylab.plot(x, t)
    117 
    118     plot2.set_ylim([-1,10])
     117    plot2.set_ylim([-1,20])
    119118    pylab.xlabel('Position')
    120119    pylab.ylabel('Height')
    121120
    122     plot3 = pylab.subplot(313)
    123 
     121    plot3 = pylab.subplot(513)
    124122    vplot, = pylab.plot(x, v)
    125 
    126     plot3.set_ylim([-15,15])
     123    plot3.set_ylim([-30,30])
    127124    pylab.xlabel('Position')
    128125    pylab.ylabel('Velocity')
     126
     127    plot4 = pylab.subplot(514)
     128    splot, = pylab.plot(x, s)
     129    plot4.set_ylim([-1,2])
     130    pylab.xlabel('Position')
     131    pylab.ylabel('State')
     132
     133    plot5 = pylab.subplot(515)
     134    Mplot, = pylab.plot(x, m)
     135    mplot, = pylab.plot(x, M)
     136    plot5.set_ylim([-1,45000])
     137    pylab.xlabel('Position')
     138    pylab.ylabel('Mass')
     139
    129140
    130141    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     
    134145        t = domain.get_quantity('top', 'centroids')
    135146        v = domain.get_quantity('velocity', 'centroids')
    136 
     147        s = domain.state
     148        m = domain.get_mass()
     149        M = m.sum() * numpy.ones_like(x)
     150       
    137151        zplot.set_ydata(z)
    138152        ztplot.set_ydata(z+t)
     
    141155        tplot.set_ydata(t)
    142156        vplot.set_ydata(v)
    143 
     157        splot.set_ydata(s)
     158        mplot.set_ydata(m)
     159        Mplot.set_ydata(M)
     160       
    144161        pylab.draw()
    145162
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

    r8263 r8265  
    264264            #Pass control on to outer loop for more specific actions
    265265            yield(t)
     266
     267
     268    def get_mass(self):
     269        """Returns array of mass in cells
     270        """
     271
     272        # Equivalent area A is defined by A = (\rho/\rho_0) A_{pipe}
     273        # A_{pipe} = width * top
     274        # A = width * height
     275        # \rho_0 = 1
     276        # Therefore \rho = height/top for pressurised states
     277        # For free surface \rho = 1
     278
     279        # Coarse approximation of mass M in a cell, assume rho is constant
     280        # M = \Delta x * \rho * A
     281        #   = \Delta_x * A * height/top (pressurised)
     282        #   = \Delta_x * A (un-pressurised)
     283
     284        area      = self.quantities['area']       
     285        height    = self.quantities['height']       
     286        top       = self.quantities['top']
     287       
     288        #Arrays   
     289        a   = area.centroid_values       
     290        h   = height.centroid_values
     291        t   = top.centroid_values
     292       
     293        M = numpy.where(self.state == 0, self.areas * a, self.areas * a * h/t)
     294
     295        return M
    266296
    267297# Auxillary methods
Note: See TracChangeset for help on using the changeset viewer.