Changeset 3362

Timestamp:

Jul 18, 2006, 9:39:13 PM (18 years ago)

Author:

jakeman

Message:

New files dam2.py and dry_dam.py further test numerical solution
against analytical solution of
Stoker 57. Limiting is still not working correctly.

Location:

development/pyvolution-1d

Files:

• dam.py (modified) (3 diffs)
• dam2.py (added)
• dry_dam.py (added)
• quantity.py (modified) (4 diffs)
• shallow_water_1d.py (modified) (5 diffs)

development/pyvolution-1d/dam.py

@@ +1 @@
+
 import os
 from math import sqrt, pi
@@ -84 +85 @@
 
 def stage(x):
-    y = zeros(len(x))
+    y = zeros(len(x),Float)
     for i in range(len(x)):
         if x[i]<=1000.0:
@@ -123 +124 @@
         StageQ = domain.quantities['stage'].vertex_values
         y = analytical_sol(X.flat,domain.time)
-        #linePlot(plot1,X,StageQ,X,y)
+        linePlot(plot1,X,StageQ,X,y)
     #raw_input('press_return')
     #pass

development/pyvolution-1d/quantity.py

@@ -513 +513 @@
         """
 
+        Z = self.gradients
+        #print "gradients 1",Z
+        
         self.compute_gradients()
+        #print "gradients 2", Z
 
         G = self.gradients
@@ -527 +531 @@
             #vertex coordinates
             x0, x1 = V[k,:]
+
             #Extrapolate
             Qv[k,0] = Qc[k] + G[k]*(x0-x)
             Qv[k,1] = Qc[k] + G[k]*(x1-x)
 
-    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
 
@@ -542 +547 @@
         postcondition:
         vertex values are updated
-        """
+        
         print "in Q.limit"
         from Numeric import zeros, Float
@@ -591 +596 @@
             for i in range(2):
                 qv[k,i] = qc[k] + phi*dq[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
+        """
+
+        from Numeric import zeros, Float
+
+        N = self.domain.number_of_elements
+        beta = self.beta
+
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        #Find min and max of this and neighbour's centroid values
+        qmax = self.qmax
+        qmin = self.qmin
+
+        for k in range(N):
+            qmax[k] = qmin[k] = qc[k]
+            for i in range(2):
+                n = self.domain.neighbours[k,i]
+                if n >= 0:
+                    qn = qc[n] #Neighbour's centroid value
+
+                    qmin[k] = min(qmin[k], qn)
+                    qmax[k] = max(qmax[k], qn)
+
+
+        #Diffences between centroids and maxima/minima
+        dqmax = qmax - qc
+        dqmin = qmin - qc
+
+        #Deltas between vertex and centroid values
+        dq = zeros(qv.shape, Float)
+        for i in range(2):
+            dq[:,i] = qv[:,i] - qc
+
+        #Phi limiter
+        for k in range(N):
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(2):
+                r = 1.0
+                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
+                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
+
+                phi = min( min(r*beta, 1), phi )
+
+            #Then update using phi limiter
+            for i in range(2):
+                qv[k,i] = qc[k] + phi*dq[k,i]
+
+        
 
 def newLinePlot(title='Simple Plot'):

development/pyvolution-1d/shallow_water_1d.py

@@ -1 @@
-"""Class Domain -
+u"""Class Domain -
 1D interval domains for finite-volume computations of
 the shallow water wave equation.
@@ -165 +165 @@
         #self.check_integrity()
 
-        msg = 'Parameter beta_h must be in the interval [0, 1)'
-        assert 0 <= self.beta_h < 1.0, msg
-        msg = 'Parameter beta_w must be in the interval [0, 1)'
-        assert 0 <= self.beta_w < 1.0, msg
+        #msg = 'Parameter beta_h must be in the interval [0, 1)'
+        #assert 0 <= self.beta_h < 1.0, msg
+        #msg = 'Parameter beta_w must be in the interval [0, 1)'
+        #assert 0 <= self.beta_w < 1.0, msg
 
 
@@ -605 +605 @@
 
     #Remove very thin layers of water
-    #protect_against_infinitesimal_and_negative_heights(domain)
+    protect_against_infinitesimal_and_negative_heights(domain)
 
     #Extrapolate all conserved quantities
@@ -634 +634 @@
 
     #Take bed elevation into account when water heights are small
-    #balance_deep_and_shallow(domain)
+    balance_deep_and_shallow(domain)
 
     #Compute edge values by interpolation
@@ -662 +662 @@
 
             #Control momentum
-#            xmomc[k] = ymomc[k] = 0.0
+            #xmomc[k] = ymomc[k] = 0.0
             xmomc[k] = 0.0