Changeset 8265
- Timestamp:
- Dec 6, 2011, 11:09:34 AM (12 years ago)
- 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 100 100 v = domain.get_quantity('velocity', 'centroids') 101 101 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) 105 107 zplot, = pylab.plot(x, z) 106 108 wplot, = pylab.plot(x, w) 107 ztplot, = pylab.plot(x, z+t) 108 109 ztplot, = pylab.plot(x, z+t) 109 110 plot1.set_ylim([-1,50]) 110 111 pylab.xlabel('Position') 111 112 pylab.ylabel('Stage') 112 113 113 plot2 = pylab.subplot(312) 114 114 plot2 = pylab.subplot(512) 115 115 hplot, = pylab.plot(x, h) 116 116 tplot, = pylab.plot(x, t) 117 118 plot2.set_ylim([-1,10]) 117 plot2.set_ylim([-1,20]) 119 118 pylab.xlabel('Position') 120 119 pylab.ylabel('Height') 121 120 122 plot3 = pylab.subplot(313) 123 121 plot3 = pylab.subplot(513) 124 122 vplot, = pylab.plot(x, v) 125 126 plot3.set_ylim([-15,15]) 123 plot3.set_ylim([-30,30]) 127 124 pylab.xlabel('Position') 128 125 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 129 140 130 141 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): … … 134 145 t = domain.get_quantity('top', 'centroids') 135 146 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 137 151 zplot.set_ydata(z) 138 152 ztplot.set_ydata(z+t) … … 141 155 tplot.set_ydata(t) 142 156 vplot.set_ydata(v) 143 157 splot.set_ydata(s) 158 mplot.set_ydata(m) 159 Mplot.set_ydata(M) 160 144 161 pylab.draw() 145 162 -
trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py
r8263 r8265 264 264 #Pass control on to outer loop for more specific actions 265 265 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 266 296 267 297 # Auxillary methods
Note: See TracChangeset
for help on using the changeset viewer.