Changeset 3510

Timestamp:

Aug 21, 2006, 10:03:42 AM (18 years ago)

Author:

jakeman

Message:

Files now contain inmplementation of method in which advection terms and
pressure terms are split. problems are present.

Location:

development/pyvolution-1d

Files:

• domain.py (modified) (1 diff)
• quantity.py (modified) (5 diffs)
• shallow_water_1d.py (modified) (9 diffs)

development/pyvolution-1d/domain.py

@@ -22 +22 @@
 
         self.beta = 1.0
+        #self.limiter = "minmod_kurganov"
+        #self.limiter_type = "steven"
+        self.fluxfunc = "unsplit"
+        
         #Store Points
         self.coordinates = array(coordinates)

development/pyvolution-1d/quantity.py

@@ -51 +51 @@
         #self.edge_values = zeros((N, 2), Float)
         #edge values are values of the ends of each interval
-        
-        #does oe dimension need edge values???
-
+      
         #Intialise centroid values
         self.interpolate()
@@ -414 +412 @@
         self.centroid_values += timestep*self.explicit_update
 
+        #SECOND ORDER RUNGE_KUTTA
+        #Re-evaluate source term with updated centroid values
+        #and update conserved quantities accordingly
+        """
+        self.explicit_update = zeros(N, Float)
+        self.domain.compute_forcing_terms()
+        self.centroid_values += 0.5*timestep*(explicit_update
+        """
+        
 ##         #Semi implicit updates
 ##         denominator = ones(N, Float)-timestep*self.semi_implicit_update
@@ -679 +686 @@
     """
 
-    def limit(self):
-        """Limit slopes for each volume to eliminate artificial variance
+    """def limit(self):
+        Limit slopes for each volume to eliminate artificial variance
         introduced by e.g. second order extrapolator
 
@@ -690 +697 @@
         postcondition:
         vertex values are updated
-        """
-
+        
         from Numeric import zeros, Float
 
         N = self.domain.number_of_elements
-        #beta = self.beta
-        beta = 0.8
+        beta = self.beta
+        #beta = 0.8
 
         qc = self.centroid_values
@@ -741 +747 @@
                 qv[k,i] = qc[k] + phi*dq[k,i]
 
-        
-
+    """
+   
+    
+    def limit(self):
+        """Limit slopes for each volume to eliminate artificial variance
+         introduced by e.g. second order extrapolator
+         
+        This is an unsophisticated limiter as it does not take into
+        account dependencies among quantities.
+
+        precondition:
+        vertex values are estimated from gradient
+        postcondition:
+        vertex values are updated
+
+        """
+        import sys
+        from Numeric import zeros, Float
+        from util import minmod, minmod_kurganov, maxmod, vanleer
+
+        N = self.domain.number_of_elements
+        limiter = self.domain.limiter
+        limiter_type = self.domain.limiter_type
+            
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        #Find min and max of this and neighbour's centroid values
+        beta_p = zeros(N,Float)
+        beta_m = zeros(N,Float)
+        beta_x = zeros(N,Float)
+        C = self.domain.centroids
+        X = self.domain.vertices
+
+        for k in range(N):
+        
+            n0 = self.domain.neighbours[k,0]
+            n1 = self.domain.neighbours[k,1]
+            if limiter_type == "slope":
+                if ( n0 >= 0) & (n1 >= 0):
+                    #SLOPE DERIVATIVE LIMIT
+                    beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
+                    beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
+                    beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
+            elif limiter_type == "flux":
+                if (n0 >= 0) & (n1 >= 0):
+                    # Check denominator not zero
+                    if (qc[k+1]-qc[k]) == 0.0:
+                        beta_p[k] = float(sys.maxint)
+                        beta_m[k] = float(sys.maxint)
+                    else:
+                        #STEVE LIMIT
+                        beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
+                        beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
+
+        #Deltas between vertex and centroid values
+        dq = zeros(qv.shape, Float)
+        for i in range(2):
+            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
+            
+        #Phi limiter
+        for k in range(N):
+            if limiter_type == "slope":
+                if limiter == "minmod":
+                    phi = minmod(beta_p[k],beta_m[k])
+
+                elif limiter == "minmod_kurganov":
+                    # Also known as monotonized central difference limiter
+                    # if theta = 2.0
+                    theta = 2.0 
+                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
+                elif limiter == "superbee":
+                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
+                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
+                    phi = maxmod(slope1,slope2)
+                    #phi = max(0.0,min(1.0,2.0*beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0*beta_p[k]),min(2.0,beta_p[k]))-1.0
+
+                elif limiter == "vanleer":
+                    phi = vanleer(beta_p[k],beta_m[k])
+
+
+                elif limiter == "muscl":
+
+                    #MINMOD
+                    #phi = minmod(1.0,beta_m[k])
+                    #SUPERBEE
+                    #phi = max(0.0,min(2.0*beta_m[k],1.0),min(beta_m[k],2.0))
+                    #VAN LEER
+                    phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))
+                    
+                for i in range(2):
+                    qv[k,i] = qc[k] + phi*dq[k,i]
+                
+            elif limiter_type == "flux":
+                phi = 0.0
+                if limiter == "flux_minmod":
+                    #FLUX MINMOD
+                    phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
+                elif limiter == "flux_superbee":
+                    #FLUX SUPERBEE
+                    phi = max(0.0,min(1.0,2.0*beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0*beta_p[k]),min(2.0,beta_p[k]))-1.0
+                elif limiter == "flux_muscl":
+                    #FLUX MUSCL
+                    phi = max(0.0,min(2.0,2.0*beta_m[k],2.0*beta_p[k],0.5*(beta_m[k]+beta_p[k])))
+                elif limiter == "flux_vanleer":
+                    #FLUX VAN LEER
+                    phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))+(beta_p[k]+abs(beta_p[k]))/(1.0+abs(beta_p[k]))-1.0
+
+                #Then update using phi limiter
+                n = self.domain.neighbours[k,1]
+                if n>=0:
+                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
+                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
+                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
+                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
+                else:
+                    qv[k,i] = qc[k]
+                    
+            
 def newLinePlot(title='Simple Plot'):
     import Gnuplot

development/pyvolution-1d/shallow_water_1d.py

@@ -64 +64 @@
         #forcing terms not included in 1d domain ?WHy?
         self.forcing_terms.append(gravity)
-        self.forcing_terms.append(manning_friction)
+        #self.forcing_terms.append(manning_friction)
         #print "\nI have Removed forcing terms line 64 1dsw"
 
@@ -283 +283 @@
 
 
-# Flux computation
-#def flux_function(normal, ql, qr, zl, zr):
-def flux_function(normal, ql, qr, zl, zr):
+def flux_function_unsplit(normal, ql, qr, zl, zr):
     """Compute fluxes between volumes for the shallow water wave equation
     cast in terms of w = h+z using the 'central scheme' as described in
@@ -376 +374 @@
 
     return edgeflux, max_speed
+
+def flux_function_split(normal, ql, qr, zl, zr):
+    from config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #soundspeed_left  = sqrt(g*h_left)
+    #soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    s_max = max(u_left, u_right, 0)
+    
+    #Minimal wave speed
+    #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    s_min = min(u_left, u_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left*h_left])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right*h_right])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left])# + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right])# + 0.5*g*h_right*h_right])
+    
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed    
+
 
 def compute_fluxes(domain):
@@ -471 +549 @@
             #print 'qr',qr
             
-            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
+            if domain.fluxfunc == "unsplit":
+                edgeflux, max_speed = flux_function_unsplit(normal, ql, qr, zl, zr)
+            elif domain.fluxfunc == "split":
+                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
             #print 'edgeflux', edgeflux
 
@@ -480 +562 @@
             try:
                 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
-                #timestep = 0.01
-            
                 timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
-                if (timestep < 1e-6) & (enter == True):
-                    #print "domain.order", domain.order
-                    #domain.write_time()
-                    print "cell number", k
-                    print "timestep", timestep
-                    print 'max_speed', max_speed
-                    s = domain.quantities['stage']
-                    s = s.centroid_values
-                    xm = domain.quantities['xmomentum']
-                    xm = xm.centroid_values
-                    print 'h', s[k]
-                    print 'xm', xm[k]
-                    print 'u', xm[k]/s[k]
-                    enter = False
-            
             except ZeroDivisionError:
                 pass
@@ -504 +569 @@
         #quantities get updated
         flux /= domain.areas[k]
-        # ADD ABOVE LINE AGAIN
+
         Stage.explicit_update[k] = flux[0]
         Xmom.explicit_update[k] = flux[1]
@@ -608 +673 @@
 
     #Remove very thin layers of water
-    ##protect_against_infinitesimal_and_negative_heights(domain)
+    protect_against_infinitesimal_and_negative_heights(domain)
 
     #Extrapolate all conserved quantities
@@ -1001 +1066 @@
     #h = Stage.edge_values - Elevation.edge_values
     h = Stage.vertex_values - Elevation.vertex_values
-    v = Elevation.vertex_values
+    b = Elevation.vertex_values
+    w = Stage.vertex_values
 
     x = domain.get_vertex_coordinates()
@@ -1014 +1080 @@
         x0, x1 = x[k,:]
         #z0, z1, z2 = v[k,:]
-        z0, z1 = v[k,:]
+        b0, b1 = b[k,:]
+
+        w0, w1 = w[k,:]
+        wx = gradient(x0, x1, w0, w1)
 
         #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
-        zx = gradient(x0, x1, z0, z1)
+        bx = gradient(x0, x1, b0, b1)
         
         #Update momentum
-        xmom[k] += -g*zx*avg_h
+        if domain.fluxfunc == "unsplit":
+            xmom[k] += -g*bx*avg_h
+        elif domain.fluxfunc == "split":
+            xmom[k] += -g*wx*avg_h
 #        ymom[k] += -g*zy*avg_h