Changeset 1586

Timestamp:

Jul 7, 2005, 12:04:36 PM (19 years ago)

Author:

steve

Message:

 

File:

• inundation/ga/storm_surge/analytical solutions/Analytical_solution_contracting_channel.py (modified) (7 diffs)

inundation/ga/storm_surge/analytical solutions/Analytical_solution_contracting_channel.py

@@ -28 +28 @@
 from shallow_water import Constant_height, Domain
 from pmesh2domain import pmesh_to_domain_instance
-from mesh_factory import contracting_channel
+from mesh_factory import contracting_channel_cross
 
 #-------
@@ -43 +43 @@
 W_downstream = 2.5
 L_1 = 5.
-L_2 = 10.
+L_2 = 11
 L_3 = Total_length - L_1 - L_2
 n = 20
@@ -49 +49 @@
 
 points, elements, boundary = \
-     contracting_channel(m, n, W_upstream, W_downstream, L_1, L_2, L_3)
+     contracting_channel_cross(m, n, W_upstream, W_downstream, L_1, L_2, L_3)
 domain = Domain(points, elements, boundary)
 
@@ -66 +66 @@
 domain.filename = 'contracting_channel_second-order'
 
+
 #----------------------------------------------------------
 # Decide which quantities are to be stored at each timestep
 domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
+#------------------------------------------------
 # This is for Visual Python
-domain.visualise = True
+#domain.visualise = True
 
 #------------------------------------------
@@ -93 +95 @@
 
 # Use the inscribed circle with safety factor of 0.9 to establish the time step
-#domain.set_to_inscribed_circle(safety_factor=0.9)
-
-### List all triangles that exist in region x > 35
-##n_triangle = []
-##for n in range(len(domain.triangles)):
-##    t = domain.triangles[n]
-##    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
-##
-##    for i in range(n_points):
-##        if inside_polygon(points(i),tri):
-##            print 'Point is within triangle with vertices '+'%s'%tri
-##            n_triangle[i] = n
-##            t = domain.triangles[n_triangle[i]]
-##            tri[i] = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+# domain.set_to_inscribed_circle(safety_factor=0.9)
 
 #----------
@@ -112 +101 @@
 import time
 t0 = time.time()
-for t in domain.evolve(yieldstep = 0.2, finaltime = 50.0):
+for t in domain.evolve(yieldstep = 5.0, finaltime = 50.0):
     domain.write_time()
 
@@ -131 +120 @@
         n_points = n_points + 1
 
-
 average_stage = average_stage/n_points
 
 #Standard Deviation
 sigma = 0.0
+max_stage = -999999.
+min_stage = 999999
 for n in range(N):
     if XY[n,0] > 35.0:
         sigma = sigma + (average_stage - stage[n])**2
+        if stage[n] > max_stage:
+            max_stage = stage[n]
+        if stage[n] < min_stage:
+            min_stage = stage[n]
 
 import math
 sigma = math.sqrt(sigma/(n_points-1))
 
-print average_stage, sigma
+print L_2, average_stage, sigma, max_stage, min_stage, n_points
 
-
-
-
-
-##average_stage = 0.
-##t_array = asarray([[0,1,2]])
-##for i in range(n_points):
-##    tri_array = asarray(tri(i))
-##    stage = domain.get_quantity('stage')[n_triangle(i)]
-##    average_triangle_stage[i] = (stage[0]+stage[1]+stage[2])/3.
-##    average_stage = average_stage + average_triangle_stage[i]
-##
-##average_stage = average_stage/n_points
-##sigma = 0.
-##for i in range(n_points):
-##    sigma = sigma + average_triangle_stage[i] - average_stage
-##
-##sigma = sqrt(sigma/n_points)
-##print sigma
-##
-
-
+domain.initialise_visualiser()
+domain.visualiser.update_all()