Changeset 989

Timestamp:

Mar 2, 2005, 5:49:23 PM (20 years ago)

Author:

steve

Message:

Corrected h limited update to ensure conservation.

Location:

inundation/ga/storm_surge

Files:

• analytical solutions/Analytical_solution_Yoon_import_mesh.py (modified) (4 diffs)
• pyvolution-1d/quantity.py (modified) (1 diff)
• pyvolution/shallow_water.py (modified) (4 diffs)
• pyvolution/shallow_water_ext.c (modified) (3 diffs)
• pyvolution/util.py (modified) (1 diff)

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

@@ -4 +4 @@
 and used by Yoon and Chou.
 
-   Copyright 2004
-   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
-   Geoscience Australia
+   Copyright 2005
+   Christopher Zoppou, Stephen Roberts, Ole Nielsen
+   ANU, Geoscience Australia
    
 """
 
-######################
+#---------------
 # Module imports 
 import sys
 from os import sep
 sys.path.append('..'+sep+'pyvolution')
-
 from shallow_water import Domain, Dirichlet_boundary
 from math import sqrt, cos, sin, pi
 from mesh_factory import strang_mesh
+from util import inside_polygon
+from Numeric import asarray
+from least_squares import Interpolation
 
-
-######################
-# Domain
-# Strang_domain will search through the file and test to see if there are
-# two or three entries. Two entries are for points and three for triangles.
-
+#-------------------------------
+# Set up the domain of triangles
+# Strang_domain will search
+# through the file and test to
+# see if there are two or three
+# entries. Two entries are for
+# points and three for triangles.
 points, elements = strang_mesh('yoon_circle.pt')
 domain = Domain(points, elements) 
 
+#----------------
+# Order of scheme
 domain.default_order = 2
+
 domain.smooth = True
 
-
+#-------------------------------------
 # Provide file name for storing output
-domain.store = True
+domain.store = False
 domain.format = 'sww'
-domain.filename = 'yoon_strang_second_order'
-
+domain.filename = 'yoon_mesh_second_order.2'
 print 'Number of triangles = ', len(domain)
 
-#Reduction operation for get_vertex_values
+#----------------------------------------------------------
+# Decide which quantities are to be stored at each timestep
+domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+#------------------------------------------
+# Reduction operation for get_vertex_values
 from util import mean
-domain.reduction = mean 
-#domain.reduction = min  #Looks better near steep slopes
+domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
 
 
-######################
-#Initial condition
-#
+#------------------
+# Initial condition
 print 'Initial condition'
 t = 0.0
@@ -58 +66 @@
 A = (L**4 - R0**4)/(L**4 + R0**4)
 omega = 2./L*sqrt(2.*g*D0)
-print 'omega = ',omega
 T = pi/omega
-print T
 
-
-#Set bed-elevation and friction(None)
+#------------------
+# Set bed elevation
 def x_slope(x,y):
     n = x.shape[0]
@@ -71 +77 @@
         z[i] = -D0*(1.-r*r/L/L)
     return z
-
 domain.set_quantity('elevation', x_slope)
 
-#Set the initial water stage
-def stage(x,y):
+#----------------------------
+# Set the initial water level
+def level(x,y):
     z = x_slope(x,y)
     n = x.shape[0]
@@ -86 +92 @@
         h[i] = z[i]
     return h   
-
-domain.set_quantity('stage', stage)
+domain.set_quantity('stage', level)
 
 
-############
-#Boundary
+#---------
+# Boundary
+print 'Boundary conditions'
 domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
+
+#---------------------------------------------
+# Find triangle that contains the point points
+# and print to file
+points = [0.,0.]
+for n in range(len(domain.triangles)):
+    t = domain.triangles[n]
+    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+    if inside_polygon(points,tri):
+        print 'Point is within triangle with vertices '+'%s'%tri
+        n_point = n
+t = domain.triangles[n_point]
+tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+filename=domain.filename
+file = open(filename,'w')
     
-
-######################
-#Evolution
+#----------
+# Evolution
 import time
 t0 = time.time()
-for t in domain.evolve(yieldstep = 10.0, finaltime = 5000):
+for t in domain.evolve(yieldstep = 20.0, finaltime = 3000 ):
     domain.write_time()
 
+    tri_array = asarray(tri)
+    t_array = asarray([[0,1,2]])
+    interp = Interpolation(tri_array,t_array,[points])
+
+    stage = domain.get_quantity('stage')[n_point]
+    xmomentum = domain.get_quantity('xmomentum')[n_point]
+    ymomentum = domain.get_quantity('ymomentum')[n_point]
+
+    interp_stage = interp.interpolate([[stage[0]],[stage[1]],[stage[2]]])
+    interp_xmomentum = interp.interpolate([[xmomentum[0]],[xmomentum[1]],[xmomentum[2]]])
+    interp_ymomentum = interp.interpolate([[ymomentum[0]],[ymomentum[1]],[ymomentum[2]]])
+    
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,interp_stage[0][0],interp_xmomentum[0][0],interp_ymomentum[0][0]) )
+
+file.close()    
 print 'That took %.2f seconds' %(time.time()-t0)
-
-    

inundation/ga/storm_surge/pyvolution-1d/quantity.py

@@ -446 +446 @@
     Q1.limit()
     
-    from matplotlib.matlab import *
+    from pylab import *
     plot(Xc,Qc)
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -59 +59 @@
 Generic_domain = Domain #Rename
 
+#Shalow water domain
 class Domain(Generic_domain):
 
@@ -435 +436 @@
             #Update optimal_timestep
             try:
-                timestep = min(timestep, domain.radii[k]/max_speed)
+                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
             except ZeroDivisionError:
                 pass
@@ -1501 +1502 @@
 
 
+
+def compute_fluxes_python(domain):
+    """Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities 
+    stage, xmomentum and ymomentum    
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes. 
+    """
+
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Bed = domain.quantities['elevation']        
+
+    #Arrays
+    stage = Stage.edge_values
+    xmom =  Xmom.edge_values
+    ymom =  Ymom.edge_values
+    bed =   Bed.edge_values    
+
+    stage_bdry = Stage.boundary_values
+    xmom_bdry =  Xmom.boundary_values
+    ymom_bdry =  Ymom.boundary_values
+    
+    flux = zeros((N,3), Float) #Work array for summing up fluxes
+
+    #Loop
+    timestep = float(sys.maxint)    
+    for k in range(N):
+
+        for i in range(3):
+            #Quantities inside volume facing neighbour i
+            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
+            zl = bed[k, i]
+
+            #Quantities at neighbour on nearest face
+            n = domain.neighbours[k,i] 
+            if n < 0:
+                m = -n-1 #Convert negative flag to index
+                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
+                zr = zl #Extend bed elevation to boundary
+            else:    
+                m = domain.neighbour_edges[k,i]
+                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
+                zr = bed[n, m]                
+
+                
+            #Outward pointing normal vector   
+            normal = domain.normals[k, 2*i:2*i+2]
+
+            #Flux computation using provided function
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+        
+            flux[k,:] = edgeflux
+
+    return flux
+    
+  
+
+
+
+
+
 ##############################################
 #Initialise module
@@ -1519 +1601 @@
     
     #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
+
 
 

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -626 +626 @@
       //timestep = min(timestep, domain.radii[k]/max_speed)
       if (max_speed > epsilon) {
-        timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+        timestep = min(timestep, 0.5*((double *) radii -> data)[k]/max_speed);
       }   
     } // end for i
@@ -767 +767 @@
   hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
 
-      
+  
+
   //Find min and max of this and neighbour's centroid values
   hmin = malloc(N * sizeof(double));
@@ -779 +780 @@
     for (i=0; i<3; i++) {
       n = ((long*) neighbours -> data)[k3+i];
+      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k];  
       if (n >= 0) {
         hn = ((double*) hc -> data)[n]; //Neighbour's centroid value

inundation/ga/storm_surge/pyvolution/util.py

@@ -1016 +1016 @@
        Polygon_function(polygons)
        
-    where polygons is a dictionary of the form
-    
-      {P0: v0, P1: v1, ...}
+    where polygons is a list of tuples of the form
+    
+      [ (P0, v0), (P1, v1), ...]
       
       with Pi being lists of vertices defining polygons and vi either