Changeset 9230

Timestamp:

Jul 1, 2014, 11:51:49 AM (10 years ago)

Author:

steve

Message:

Updating supercritical overbump test

Location:

trunk/anuga_core/validation_tests/analytical_exact

Files:

• subcritical_over_bump/produce_results.py (modified) (1 diff)
• supercritical_over_bump/analytical_supercritical.py (moved) (moved from trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/analytical_subcritical.py) (3 diffs)
• supercritical_over_bump/numerical_supercritical.py (moved) (moved from trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/numerical_subcritical.py) (5 diffs)
• supercritical_over_bump/plot_results.py (modified) (3 diffs)

trunk/anuga_core/validation_tests/analytical_exact/subcritical_over_bump/produce_results.py

@@ -4 +4 @@
 args = anuga.get_args()
 
-produce_report('numerical_subcritical.py', args=args)
+produce_report('numerical_supercritical.py', args=args)
 
 

trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/analytical_supercritical.py

@@ -1 +1 @@
 """
-Subcritical flow over a bump.
-Ref1: Houghton & Kasahara, Nonlinear shallow fluid flow over an isolated ridge.
-Comm. Pure and Applied Math. DOI:10.1002/cpa.3160210103
+Supercritical flow over a bump.
 
-Ref2: Delestre et al, 2012, SWASHES: a compilation of shallow water
-analytic solutions..., Int J Numer Meth Fluids, DOI:10.1002/fld.3741
 
-Sudi Mungkasi, ANU 2012
+Steve Roberts, ANU 2014
 """
-from numpy import zeros, linspace
+from numpy import zeros
 from scipy.optimize import fsolve
-from pylab import plot, ylim, show
 from anuga import g
 
-qA = 4.42  # This is the imposed momentum
-hx = 2.0   # This is the water height downstream
+
+qA = 10.0  # This is the imposed momentum
+hx = 0.5   # This is the water height downstream
 
 def analytic_sol(x):    
+
     def elevation(x):
         z_b = zeros(len(x))
@@ -28 +25 @@
     z = elevation(x)
 
+ 
     def find_hL(h): #to find the water height at every spatial point
         return h**3 + (zb-qA**2/(2*g*hx**2)-hx)*h**2 + qA**2/(2*g)
@@ -33 +31 @@
     for i in range(len(x)):
         zb = z[i]
-        h[i] = fsolve(find_hL, 2.0)
+        h[i] = fsolve(find_hL, 0.5)
+
+
     return h,z
 
-##N = 401
-##L = 25.
-##x = linspace(0.0,L,N)
-##h,z = analytic_sol(x)
-##plot(x,h+z, x,z)
-##ylim([-0.1, 2.1])
-##show()
+

trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/numerical_supercritical.py

@@ -1 +1 @@
 """
 Simple water flow example using ANUGA
-Subcritical flow over a bump.
+Supercritical flow over a bump.
 """
 
@@ -25 +25 @@
 #output_dir = 'subcritical_'+time
 output_dir = '.'
-output_file = 'subcritical'
+output_file = 'supercritical'
 
 args = anuga.get_args()
@@ -35 +35 @@
 # Setup domain
 #------------------------------------------------------------------------------
-dx = 0.1
+dx = 0.05
 dy = dx
 L = 25.
@@ -80 +80 @@
 # Setup boundary conditions
 #------------------------------------------------------------------------------
-from math import sin, pi, exp
 Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
 #Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
-Bd = anuga.Dirichlet_boundary([2., 4.42, 0.]) # Constant boundary values
+Bd = anuga.Dirichlet_boundary([0.5, 10.0, 0.]) # Constant boundary values
 
 # Associate boundary tags with boundary objects
@@ -105 +104 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-for t in domain.evolve(yieldstep = 1.0, finaltime = 100.):
+for t in domain.evolve(yieldstep = 1.0, finaltime = 10.):
     #print domain.timestepping_statistics(track_speeds=True)
     if myid == 0 and verbose: print domain.timestepping_statistics()

trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/plot_results.py

@@ -3 +3 @@
 
 """
+
 import anuga.utilities.plot_utils as util
 from matplotlib import pyplot as pyplot
-from analytical_subcritical import *
+from analytical_supercritical import analytic_sol
 from numpy import ones, arange
 
-p_st = util.get_output('subcritical.sww')
+import anuga
+parser = anuga.create_standard_parser()
+
+parser.add_argument('-tid', type=int, default=-1, help ='timestep id')
+args = parser.parse_args()
+tid = args.tid
+verbose = args.verbose
+
+
+if verbose: print 'Read in swwfile'
+p_st = util.get_output('supercritical.sww')
 p2_st=util.get_centroids(p_st)
 
-v = p2_st.y[10]
+#v = p2_st.y[10]
 #v2=(p2_st.y==v)
-v#2=(p2_st.y>-1.0)
+#v2=(p2_st.y>-1.0)
 v2 = arange(len(p2_st.y))
 
+if verbose: print 'Calculate analytical solution'
 h,z = analytic_sol(p2_st.x[v2])
+qexact = 10
+
+
+
 
 #Plot the stages##############################################################
+if verbose: print 'Create Stage plot'
 pyplot.clf()
-pyplot.plot(p2_st.x[v2], p2_st.stage[-1,v2], 'b.-', label='numerical stage') # 0*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[tid,v2], 'b.-', label='numerical stage') # 0*T/6
 pyplot.plot(p2_st.x[v2], h+z,'r-', label='analytical stage')
 pyplot.plot(p2_st.x[v2], z,'k-', label='bed elevation')
-pyplot.title('Stage at time = %s secs'% p2_st.time[-1])
+pyplot.title('Stage at time = %s secs'% p2_st.time[tid])
 ##pyplot.ylim(-5.0,5.0)
 pyplot.legend(loc='best')
@@ -32 +49 @@
 
 #Plot the momentums##########################################################
+if verbose: print 'Create Momentum plot'
 pyplot.clf()
-pyplot.plot(p2_st.x[v2], p2_st.xmom[-1,v2], 'b.-', label='numerical') # 0*T/6
-pyplot.plot(p2_st.x[v2], 4.42*ones(len(p2_st.x[v2])),'r-', label='analytical')
-pyplot.title('Xmomentum at time = %s secs'% p2_st.time[-1])
+pyplot.plot(p2_st.x[v2], p2_st.xmom[tid,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], qexact*ones(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xmomentum at time = %s secs'% p2_st.time[tid])
 pyplot.legend(loc='best')
 ##pyplot.ylim([1.52,1.54])
@@ -45 +63 @@
 
 #Plot the velocities#########################################################
+if verbose: print 'Create Velocity plot'
 pyplot.clf()
-pyplot.plot(p2_st.x[v2], p2_st.xvel[-1,v2], 'b.-', label='numerical') # 0*T/6
-pyplot.plot(p2_st.x[v2], 4.42/h,'r-', label='analytical')
-pyplot.title('Xvelocity at time = %s secs'% p2_st.time[-1])
+pyplot.plot(p2_st.x[v2], p2_st.xvel[tid,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], qexact/h,'r-', label='analytical')
+pyplot.title('Xvelocity at time = %s secs'% p2_st.time[tid])
 pyplot.legend(loc='best')
 pyplot.xlabel('Xposition')