Changeset 335

Timestamp:

Sep 21, 2004, 5:43:03 PM (20 years ago)

Author:

steve

Message:

Limit and gradient aded to quantity

Location:

inundation/ga/storm_surge

Files:

• pyvolution-1d/domain.py (modified) (1 diff)
• pyvolution-1d/quantity.py (modified) (6 diffs)
• pyvolution-1d/test_quantity.py (modified) (5 diffs)
• pyvolution/netherlands.py (modified) (2 diffs)

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

@@ -17 +17 @@
 
         from Numeric import array, zeros, Float, Int
-        
+
+        self.beta = 1.0
         #Store Points
         self.coordinates = array(coordinates)

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

@@ -200 +200 @@
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
+
+        self.gradients = zeros(N, Float)
+        self.qmax = zeros(self.centroid_values.shape, Float)
+        self.qmin = zeros(self.centroid_values.shape, Float)              
                 
 
@@ -226 +230 @@
 
     def compute_gradients(self):
-        """Compute gradients of triangle surfaces defined by centroids of
+        """Compute gradients of piecewise linear function defined by centroids of
         neighbouring volumes.
-        If one face is on the boundary, use own centroid as neighbour centroid.
-        If two or more are on the boundary, fall back to first order scheme.
-        
-        Also return minimum and maximum of conserved quantities 
         """
 
         
         from Numeric import array, zeros, Float
-        from util import gradient
 
         N = self.centroid_values.shape[0]
 
-        a = zeros(N, Float)
+
+        G = self.gradients
+        Q = self.centroid_values
+        X = self.domain.centroids
         
         for k in range(N):
@@ -246 +248 @@
             # first and last elements have boundaries
 
-            if k == 1:
+            if k == 0:
 
                 #Get data
@@ -252 +254 @@
                 k1 = k+1
                 
-                q0 = self.centroid_values[k0]
-                q1 = self.centroid_values[k1]
-
-                x0 = self.domain.centroids[k0] #V0 centroid
-                x1 = self.domain.centroids[k1] #V1 centroid        
+                q0 = Q[k0]
+                q1 = Q[k1]
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k1] #V1 centroid        
 
                 #Gradient
-                a[k] = (q1 - q0)/(x1 - x0)
+                G[k] = (q1 - q0)/(x1 - x0)
 
             elif k == N-1:
@@ -267 +269 @@
                 k1 = k-1
                 
-                q0 = self.centroid_values[k0]
-                q1 = self.centroid_values[k1]
-
-                x0 = self.domain.centroids[k0] #V0 centroid
-                x1 = self.domain.centroids[k1] #V1 centroid        
+                q0 = Q[k0]
+                q1 = Q[k1]
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k1] #V1 centroid        
 
                 #Gradient
-                a[k] = (q1 - q0)/(x1 - x0)
+                G[k] = (q1 - q0)/(x1 - x0)
                 
-            else
+            else:
                 #Interior Volume (2 neighbours)
 
                 #Get data
-                k0 = self.domain.neighbours[k,0]
-                k2 = self.domain.neighbours[k,1]
-
-                q0 = self.centroid_values[k0]
-                q1 = self.centroid_values[k]
-                q2 = self.centroid_values[k2]                    
-
-                x0 = self.domain.centroids[k0] #V0 centroid
-                x1 = self.domain.centroids[k1] #V1 centroid (Self)
-                x2 = self.domain.centroids[k2] #V2 centroid
+                k0 = k-1
+                k2 = k+1
+
+                q0 = Q[k0]
+                q1 = Q[k]
+                q2 = Q[k2]                    
+
+                x0 = X[k0] #V0 centroid
+                x1 = X[k] #V1 centroid (Self)
+                x2 = X[k2] #V2 centroid
 
                 #Gradient
-                a[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
-        
-        return a
-        
-
-##     def limit(self):
-##         #Call correct module function
-##         #(either from this module or C-extension)
-##         limit(self)
-        
-        
-##     def extrapolate_first_order(self):
-##         """Extrapolate conserved quantities from centroid to
-##         vertices for each volume using
-##         first order scheme.
-##         """
-        
-##         qc = self.centroid_values
-##         qv = self.vertex_values
-
-##         for i in range(3):
-##             qv[:,i] = qc
-            
-            
-##     def extrapolate_second_order(self):
-##         """Extrapolate conserved quantities from centroid to
-##         vertices for each volume using
-##         second order scheme.
-##         """
-        
-##         a, b = self.compute_gradients()
-
-##         V = self.domain.get_vertex_coordinates()
-##         qc = self.centroid_values
-##         qv = self.vertex_values
-        
-##         #Check each triangle
-##         for k in range(self.domain.number_of_elements):
-##             #Centroid coordinates            
-##             x, y = self.domain.centroids[k]
-
-##             #vertex coordinates
-##             x0, y0, x1, y1, x2, y2 = V[k,:]
-
-##             #Extrapolate
-##             qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
-##             qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-##             qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
-
-            
-
-## def limit(quantity):        
-##     """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 = quantity.domain.number_of_elements
-    
-##     beta = quantity.domain.beta
-        
-##     qc = quantity.centroid_values
-##     qv = quantity.vertex_values
-        
-##     #Find min and max of this and neighbour's centroid values
-##     qmax = zeros(qc.shape, Float)
-##     qmin = zeros(qc.shape, Float)        
-    
-##     for k in range(N):
-##         qmax[k] = qmin[k] = qc[k]
-##         for i in range(3):
-##             n = quantity.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(3):
-##         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(3):
-##             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(3):            
-##             qv[k,i] = qc[k] + phi*dq[k,i]
+                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+        
+        return
+        
+    def extrapolate_first_order(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using
+        first order scheme.
+        """
+        
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        for i in range(2):
+            qv[:,i] = qc
+            
+            
+    def extrapolate_second_order(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using
+        second order scheme.
+        """
+        
+        self.compute_gradients()
+
+        G = self.gradients
+        V = self.domain.vertices
+        Qc = self.centroid_values
+        Qv = self.vertex_values
+        
+        #Check each triangle
+        for k in range(self.domain.number_of_elements):
+            #Centroid coordinates            
+            x = self.domain.centroids[k]
+
+            #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
+        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.domain.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 [-1,1]:
+                n = k+i
+                if (n >= 0) & (n <= N-1):
+                    qn = qc[n] #Neighbour's centroid value
+
+                    qmin[k] = min(qmin[k], qn)
+                    qmax[k] = max(qmax[k], qn)
+
+        #Phi limiter    
+        for k in range(N):
+
+            #Diffences between centroids and maxima/minima
+            dqmax = qmax[k] - qc[k] 
+            dqmin = qmin[k] - qc[k]
+
+            #Deltas between vertex and centroid values
+            dq = [0.0, 0.0]
+            for i in range(2):
+                dq[i] = qv[k,i] - qc[k]            
+
+            #Find the gradient limiter (phi) across vertices
+            phi = 1.0
+            for i in range(2):
+                r = 1.0
+                if (dq[i] > 0): r = dqmax/dq[i]
+                if (dq[i] < 0): r = dqmin/dq[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[i]
 
 
@@ -442 +437 @@
     print Q1.vertex_values
     print Q1.centroid_values
-
-
-
-
+    Xc = Q1.domain.vertices
+    Qc = Q1.vertex_values
+    print Xc
+    print Qc
+
+    Qc[1,0] = 3
+    from matplotlib.matlab import *
+    plot(Xc,Qc)
+
+    
+    
+
+    
+
+
+
+

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

@@ -7 +7 @@
 from quantity import *
 #from config import epsilon
-from Numeric import allclose, array
+from Numeric import allclose, array, zeros, Float
 
         
@@ -97 +97 @@
         assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
 
-        #Test exceptions
+        #Test exceptionsdatamining.anu.edu.au/svn/
         try:
             quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
@@ -112 +112 @@
             raise 'should have raised AssertionError'
 
-
-
     def test_set_values_const(self):
         quantity = Quantity(self.domain5)
@@ -144 +142 @@
 
 
-##     def test_gradient(self):
-##         quantity = Conserved_quantity(self.mesh4)
-
-##         #Set up for a gradient of (3,0) at mid triangle 
-##         quantity.set_values([2.0, 4.0, 8.0, 2.0],
-##                             location = 'centroids')
-        
-##         #Gradients
-##         a, b = quantity.compute_gradients()
-
-
-##         #gradient bewteen t0 and t1 is undefined as det==0
-##         assert a[0] == 0.0
-##         assert b[0] == 0.0
-##         #The others are OK
-##         for i in range(1,4):
-##             assert a[i] == 3.0
-##             assert b[i] == 0.0
-
-
-##         quantity.extrapolate_second_order()
-        
-##         assert allclose(quantity.vertex_values, [[2., 2.,  2.],
-##                                                  [0., 6.,  6.],
-##                                                  [6., 6., 12.],
-##                                                  [0., 0.,  6.]])
-
-
-
-##     def test_second_order_extrapolation2(self):
-##         quantity = Conserved_quantity(self.mesh4)        
-
-##         #Set up for a gradient of (3,1), f(x) = 3x+y
-##         quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
-##                             location = 'centroids')
-        
-##         #Gradients
-##         a, b = quantity.compute_gradients()
-
-##         #gradient bewteen t0 and t1 is undefined as det==0
-##         assert a[0] == 0.0
-##         assert b[0] == 0.0
-##         #The others are OK
-##      163.968025   for i in range(1,4):
-##             assert allclose(a[i], 3.0)
-##             assert allclose(b[i], 1.0)
-
-
-##         quantity.extrapolate_second_order()
-
-##         assert allclose(quantity.vertex_values[1,0], 2.0)
-##         assert allclose(quantity.vertex_values[1,1], 6.0)        
-##         assert allclose(quantity.vertex_values[1,2], 8.0)
-
-
-
-## #     def test_limiter(self):
-
-## #         initialise_consecutive_datastructure(points=6+4, elements=4)          
-        
-## #         a = Point (0.0, 0.0)
-## #         b = Point (0.0, 2.0)
-## #         c = Point (2.0, 0.0)
-## #         d = Point (0.0, 4.0)
-## #         e = Point (2.0, 2.0)
-## #         f = Point (4.0, 0.0)
-
-## #         #Set up for a gradient of (3,1), f(x) = 3x+y
-## #         v1 = Volume(b,a,c,array([0.0,0,0]))        
-## #         v2 = Volume(b,c,e,array([1.0,0,0]))
-## #         v3 = Volume(e,c,f,array([10.0,0,0]))
-## #         v4 = Volume(d,b,e,array([0.0,0,0]))
-
-## #         #Setup neighbour structure
-## #         domain = Domain([v1,v2,v3,v4])
-## #         domain.precompute()
-        
-## #         #Lets's check first order first, hey
-## #    domain.order = 1
-## #    domain.limiter = None
-## #         distribute_to_vertices_and_edges(domain)
-## #         assert allclose(v2.conserved_quantities_vertex0,
-## #                         v2.conserved_quantities_centroid)
-## #         assert allclose(v2.conserved_quantities_vertex1,
-## #                         v2.conserved_quantities_centroid)
-## #         assert allclose(v2.conserved_quantities_vertex2,
-## #                         v2.conserved_quantities_centroid)
-
-
-## #         #Gradient of fitted pwl surface
-## #         a, b = compute_gradient(v2.id)
-        
-
-## #         assert abs(a[0] - 5.0) < epsilon
-## #         assert abs(b[0]) < epsilon        
-## #         #assert qminr[0] == 0.0
-## #         #assert qmaxr[0] == 10.0
-        
-## #         #And now for the second order stuff        
-## #         # - the full second order extrapolation
-## #    domain.order = 2        
-## #         distribute_to_vertices_and_edges(domain)
-
-
-## #         qmin = qmax = v2.conserved_quantities_centroid
-
-## #         qmin = minimum(qmin, v1.conserved_quantities_centroid)
-## #         qmax = maximum(qmax, v1.conserved_quantities_centroid)
-
-## #         qmin = minimum(qmin, v3.conserved_quantities_centroid)
-## #         qmax = maximum(qmax, v3.conserved_quantities_centroid)
-
-## #         qmin = minimum(qmin, v4.conserved_quantities_centroid)
-## #         qmax = maximum(qmax, v4.conserved_quantities_centroid)                    
-## #         #assert qminr == qmin
-## #         #assert qmaxr == qmax
-
-## #         assert v2.conserved_quantities_vertex0 <= qmax
-## #         assert v2.conserved_quantities_vertex0 >= qmin
-## #         assert v2.conserved_quantities_vertex1 <= qmax
-## #         assert v2.conserved_quantities_vertex1 >= qmin
-## #         assert v2.conserved_quantities_vertex2 <= qmax
-## #         assert v2.conserved_quantities_vertex2 >= qmin                    
-
-
-## #         #Check that volume has been preserved    
-
-## #         q = v2.conserved_quantities_centroid[0]
-## #         w = (v2.conserved_quantities_vertex0[0] +
-## #              v2.conserved_quantities_vertex1[0] +
-## #              v2.conserved_quantities_vertex2[0])/3
-
-## #         assert allclose(q, w)
+    def test_gradient(self):
+        
+        quantity = Conserved_quantity(self.domain5)
+
+        #Set up for a gradient of (3,0) at mid triangle 
+        quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0],
+                            location = 'centroids')
+        
+        #Gradients
+        quantity.compute_gradients()
+        N = quantity.gradients.shape[0]
+        
+        G = quantity.gradients
+        G0 = zeros(N, Float)
+        Q = quantity.centroid_values
+        X = quantity.domain.centroids
+
+        def grad0(x0,x1,x2,q0,q1,q2):
+            return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
+
+        def grad1(x0,x1,q0,q1):
+            return  (q1 - q0)/(x1 - x0)
+
+        #Check Gradients
+        G0[0] = grad1(X[0],X[1],Q[0],Q[1])
+        for i in range(1,4):
+            G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1])
+        G0[4] = grad1(X[3],X[4],Q[3],Q[4])
+
+        assert allclose(G,G0)
+
+        quantity.extrapolate_first_order()
+        
+        assert allclose(quantity.vertex_values, [[2., 2.],
+                                                 [4., 4.],
+                                                 [8., 8.],
+                                                 [2., 2.],
+                                                 [5., 5.]])
+
+        quantity.extrapolate_second_order()
+
+        V = quantity.domain.vertices
+        for k in range(quantity.domain.number_of_elements):
+            G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5
+
+        assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]],
+                                                 [4.-G0[1], 4.+G0[1]],
+                                                 [8.-G0[2], 8.+G0[2]],
+                                                 [2.-G0[3], 2.+G0[3]],
+                                                 [5.-G0[4], 5.+G0[4]]])
+
+
+
+    def test_second_order_extrapolation2(self):
+        quantity = Conserved_quantity(self.domain5)        
+
+        #Set up for a gradient of (2), f(x) = 2x
+        quantity.set_values([1., 3., 4.5, 5.6, 7.1],
+                            location = 'centroids')
+        
+        #Gradients
+        quantity.compute_gradients()
+        G = quantity.gradients
+        
+        allclose(G, 2.0)
+
+        quantity.extrapolate_second_order()
+
+        V = quantity.domain.vertices
+        Q = quantity.vertex_values
+        
+        assert allclose(Q[:,0], 2.0*V[:,0])
+        assert allclose(Q[:,1], 2.0*V[:,1])        
+
+
+
 
 
@@ -297 +238 @@
 
 
-##     def test_second_order_extrapolator(self):
-##         quantity = Conserved_quantity(self.mesh4)                        
-
-##         #Set up for a gradient of (3,0) at mid triangle 
-##         quantity.set_values([2.0, 4.0, 8.0, 2.0],
-##                             location = 'centroids')
-        
-
-
-##         quantity.extrapolate_second_order()
-##         quantity.limit()
-
-
-##         #Assert that central triangle is limited by neighbours
-##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
-##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
-        
-##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
-##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
-        
-##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
-##         assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
-
-        
-##         #Assert that quantities are conserved
-##         from Numeric import sum
-##         for k in range(quantity.centroid_values.shape[0]):
-##             assert allclose (quantity.centroid_values[k],
-##                              sum(quantity.vertex_values[k,:])/3)
-        
-        
-
-        
-
-##     def test_limiter(self):
-##         quantity = Conserved_quantity(self.mesh4)
-
-##         #Create a deliberate overshoot (e.g. from gradient computation)
-##         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-
-
-##         #Limit
-##         quantity.limit()
-
-##         #Assert that central triangle is limited by neighbours
-##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
-##         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
-        
-##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
-##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
-        
-##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
-##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
-
-
-        
-##         #Assert that quantities are conserved
-##         from Numeric import sum
-##         for k in range(quantity.centroid_values.shape[0]):
-##             assert allclose (quantity.centroid_values[k],
-##                              sum(quantity.vertex_values[k,:])/3)
+    def test_second_order_limit_extrapolator(self):
+        quantity = Conserved_quantity(self.domain5)                        
+
+        #Set up for a gradient of (3,0) at mid triangle 
+        quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0],
+                            location = 'centroids')
+        
+
+
+        quantity.extrapolate_second_order()
+        quantity.limit()
+
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.vertex_values[0,0] >= 2.0
+        assert quantity.vertex_values[0,0] <= 4.0
+        assert quantity.vertex_values[0,1] >= 2.0
+        assert quantity.vertex_values[0,1] <= 4.0
+        
+        assert quantity.vertex_values[1,0] >= 2.0
+        assert quantity.vertex_values[1,0] <= 8.0
+        assert quantity.vertex_values[1,1] >= 2.0
+        assert quantity.vertex_values[1,1] <= 8.0
+
+        assert quantity.vertex_values[2,0] >= 2.0
+        assert quantity.vertex_values[2,0] <= 8.0
+        assert quantity.vertex_values[2,1] >= 2.0
+        assert quantity.vertex_values[2,1] <= 8.0
+
+        assert quantity.vertex_values[3,0] >= 0.0
+        assert quantity.vertex_values[3,0] <= 8.0
+        assert quantity.vertex_values[3,1] >= 0.0
+        assert quantity.vertex_values[3,1] <= 8.0
+
+        assert quantity.vertex_values[4,0] >= 0.0
+        assert quantity.vertex_values[4,0] <= 2.0
+        assert quantity.vertex_values[4,1] >= 0.0
+        assert quantity.vertex_values[4,1] <= 2.0
+        
+
+        
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/2)
+        
+        
+        
+
+    def test_limiter(self):
+        quantity = Conserved_quantity(self.domain5)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+
+
+        #Limit
+        quantity.limit()
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
+        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
+        
+        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
+        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
+        
+        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
+        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
+
+
+        
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/3)
         
 

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -92 +92 @@
 
 N = 600 #Size = 720000
-N = 220
+N = 20
 N = 55
 
@@ -172 +172 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.01, finaltime = 1.0):
+for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
     domain.write_time()