Changeset 1913

Timestamp:

Oct 14, 2005, 7:57:07 AM (19 years ago)

Author:

steve

Message:

Testing Oblique shocks in analytic solutions

File:

• inundation/analytical solutions/Analytical solution_oblique_shock_01.py (modified) (9 diffs)

inundation/analytical solutions/Analytical solution_oblique_shock_01.py

@@ -31 +31 @@
 from shallow_water import Domain, Constant_height
 from shallow_water import Transmissive_boundary, Reflective_boundary,\
-     Dirichlet_boundary
+     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary
 
 from math import sqrt, cos, sin, pi
@@ -42 +42 @@
 leny = 30.
 lenx = 40.
-n = 50
-m = 60
-n=n*4
-m=m*4
+f = 4
+n = 60*f
+m = 80*f
 
 points, elements, boundary = oblique_cross(m, n, lenx, leny)
@@ -52 +51 @@
 # Order of solver
 domain.default_order=2
+domain.beta_w = 0.8
 
 # Store output
@@ -61 +61 @@
 
 # Provide file name for storing output
-domain.filename="oblique"
+domain.filename="oblique_cross"
 
 # Visualization smoothing
 domain.smooth=True
-domain.visualise=True
+#domain.visualise=True
 
 #######################
@@ -71 +71 @@
 def x_slope(x, y):
     return 0*x
+
+def shock_hb(t):
+    return 2.5
+
+def shock_hh(t):
+    return 2.5
+
+def shock_pb(t):
+    return 10
+
+def shock_pp(t):
+    return 10
+
+
+def shock_h(x,y):
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        if x[i]<0:
+            w[i] = shock_hh(0.0)
+        else:
+            w[i] = 1.0
+    return w
+
+
+def shock_p(x,y):
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        if x[i]<0:
+            w[i] = shock_pp(0.0)
+        else:
+            w[i] = 0.0
+    return w
+
 
 domain.set_quantity('elevation', x_slope)
@@ -80 +115 @@
 R = Reflective_boundary(domain)
 T = Transmissive_boundary(domain)
-D = Dirichlet_boundary([1.82, 8.49, 0.0])
-
+D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0])
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)
 
 domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
@@ -87 +122 @@
 ######################
 #Initial condition
-h = 0.5
-domain.set_quantity('stage', Constant_height(x_slope, h) )
+
+domain.set_quantity('stage', shock_h )
+domain.set_quantity('xmomentum',shock_p )
 
 
@@ -99 +135 @@
 import time
 t0 = time.time()
-for t in domain.evolve(yieldstep = 0.1, finaltime = 5):
+for t in domain.evolve(yieldstep = 0.5, finaltime = 5.0):
     domain.write_time()
-    print domain.quantities['stage'].get_values(location='centroids',indices=[3399])
-    print domain.quantities['xmomentum'].get_values(location='centroids',indices=[3399])
-    print domain.quantities['ymomentum'].get_values(location='centroids',indices=[3399])
+    id = 3399
+    print domain.quantities['stage'].get_values(location='centroids',indices=[id])
+    print domain.quantities['xmomentum'].get_values(location='centroids',indices=[id])
+    print domain.quantities['ymomentum'].get_values(location='centroids',indices=[id])
+    id = 12719
+    print domain.quantities['stage'].get_values(location='centroids',indices=[id])
+    print domain.quantities['xmomentum'].get_values(location='centroids',indices=[id])
+    print domain.quantities['ymomentum'].get_values(location='centroids',indices=[id])
 
 print 'That took %.2f seconds' %(time.time()-t0)
@@ -114 +155 @@
 #im = ImageGrab.grab()
 #im.save("ccube.eps")
-