Changeset 7868

Timestamp:

Jun 22, 2010, 5:30:32 PM (14 years ago)

Author:

steve

Message:

Added in some more c based limiters

Location:

anuga_work/development/2010-projects/anuga_1d

Files:

• base/generic_domain.py (modified) (3 diffs)
• base/generic_mesh.py (modified) (1 diff)
• base/limiters_python.py (modified) (1 diff)
• base/quantity.py (modified) (9 diffs)
• base/quantity_ext.c (added)
• base/test_quantity_1d.py (modified) (4 diffs)
• config.py (modified) (1 diff)
• sww/run_dry_dam.py (modified) (4 diffs)
• sww/sww_domain.py (modified) (5 diffs)
• utilities/util_ext.h (modified) (1 diff)

anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py

@@ -17 +17 @@
                  coordinates,
                  boundary = None,
-                 conserved_quantities = None,
-                 evolved_quantities = None,
-                 other_quantities = None,
+                 conserved_quantities = [],
+                 evolved_quantities = [],
+                 other_quantities = [],
                  tagged_elements = None):
         """
@@ -35 +35 @@
         self.number_of_elements = N = len(self.coordinates)-1
 
-        self.beta = 1.0
-        self.set_limiter("minmod_kurganov")
         self.set_CFL(CFL)
         self.set_timestepping_method(timestepping_method)
@@ -869 +867 @@
     
 
+    def get_beta(self):
+
+        warn('limiter parameter beta associated with quantity not domain')
+
+    def set_beta(self,beta):
+        """Set the same limiter beta parameter to all evolving quantities
+        """
+
+        for name in self.evolved_quantities:
+            Q = self.quantities[name]
+            Q.set_beta(beta)
+
+
     def get_limiter(self):
-        return self.limiter
+
+        warn('limiter associated with quantity not domain')
 
     def set_limiter(self,limiter):
 
-        possible_limiters = \
-        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
-
-        if limiter in possible_limiters:
-            self.limiter = limiter
-            return
-
-        msg = '%s is an incorrect limiter type.\n'% limiter
-        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
-        raise Exception, msg
+        for name in self.evolved_quantities:
+            Q = self.quantities[name]
+            Q.set_limiter(limiter)
+        
 
         

anuga_work/development/2010-projects/anuga_1d/base/generic_mesh.py

@@ -18 +18 @@
 
 
-    points  = x_0 + (x_1 - x_0)*numpy.arange(n+1,dtype=numpy.float)
+    points  = x_0 + (x_1 - x_0)*numpy.arange(n+1,dtype=numpy.float)/n
     boundary = {(0, 0): 'left', (n-1, 1): 'right'}
 

anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py

@@ -37 +37 @@
 def vanleer(a,b):
 
-    from numpy import abs, where
+    from numpy import fabs, where
 
-    return where((abs(a) + abs(b) >= 1.0e-12), (a*abs(b)+abs(a)*b)/(abs(a)+abs(b)), 0.0)
+    return where((fabs(a) + fabs(b) >= 1.0e-12), (a*fabs(b)+fabs(a)*b)/(fabs(a)+fabs(b)), 0.0)
 
 

anuga_work/development/2010-projects/anuga_1d/base/quantity.py

@@ -70 +70 @@
         #flux calculations and forcing functions
         
-        N = domain.number_of_elements
+        self.N = domain.number_of_elements
+        N = self.N
         self.explicit_update = numpy.zeros(N, numpy.float )
         self.semi_implicit_update = numpy.zeros(N, numpy.float )
@@ -78 +79 @@
         self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
 
-        self.beta = domain.beta
+        #These are taken from domain but can be set for each quantity
+        # if so desired
+        self.beta = 0.0
+        self.limiter = 'vanleer'
 
 
@@ -604 +608 @@
         return integral
 
+    def get_beta(self,beta):
+        """Get limiting parameter
+        """
+        
+        return self.beta
+
+    def set_beta(self,beta):
+        """Set limiting parameter
+        """
+
+        #Probably should test that it is not too large
+        self.beta = beta
+
+
+    def get_limiter(self):
+        return self.limiter
+
+    def set_limiter(self,limiter):
+
+        possible_limiters = \
+        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
+
+        if limiter in possible_limiters:
+            self.limiter = limiter
+            return
+
+        msg = '%s is an incorrect limiter type.\n'% limiter
+        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
+        raise Exception, msg
 
     def update(self, timestep):
@@ -795 +828 @@
 
         #Gradient
-        G[1:-1] = xmic( self.domain.beta, d1, d2 )      
+        G[1:-1] = xmic( self.beta, d1, d2 )      
 
         
@@ -818 +851 @@
         """
 
-        if self.domain.limiter == "pyvolution":
+        if self.limiter == "pyvolution":
             self.limit_pyvolution()
 
-        elif self.domain.limiter == "minmod_steve":
+        elif self.limiter == "minmod_steve":
             self.limit_minmod()
 
@@ -985 +1018 @@
         from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
 
-        limiter = self.domain.limiter
+        limiter = self.get_limiter()
         #print limiter
         
         #print 'limit_range'
-        N = self.domain.number_of_elements
+        N = self.N
         qc = self.centroid_values
         qv = self.vertex_values
@@ -1013 +1046 @@
         phi[N-1] = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2])
 
+
         if limiter == "minmod":
             phi[1:-1] = minmod(beta_p[1:-1],beta_m[1:-1])
@@ -1023 +1057 @@
 
         elif limiter == "minmod_kurganov":
-            theta = 2.0
+            theta = self.beta
             phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
 
@@ -1034 +1068 @@
             msg = 'Unknown limiter'
             raise Exception, msg
-    
+
+
+        
         qv[:,0] = qc + phi*dx[:,0]
         qv[:,1] = qc + phi*dx[:,1]
 
 
-        #Phi limiter
-#        for k in range(N):
-#            n0 = self.domain.neighbours[k,0]
-#            n1 = self.domain.neighbours[k,1]
-#            if n0 < 0:
-#                phi = (qc[k+1] - qc[k])/(xc[k+1] - xc[k])
-#            elif n1 < 0:
-#                phi = (qc[k] - qc[k-1])/(xc[k] - xc[k-1])
-#            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
-#            #    phi = 0.0
-#            else:
-#                if limiter == "minmod":
-#                    phi = minmod(beta_p[k],beta_m[k])
-#
-#                elif limiter == "minmod_kurganov":#Change this
-#                    # 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)
-#
-#                elif limiter == "vanleer":
-#                    phi = vanleer(beta_p[k],beta_m[k])
-#
-#                elif limiter == "vanalbada":
-#                    phi = vanalbada(beta_m[k],beta_p[k])
-#
-#            for i in range(2):
-#                qv[k,i] = qc[k] + phi*dx[k,i]
+
 
     def limit_steve_slope(self):    

anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

@@ -376 +376 @@
 
         assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
-        
+
+        quantity.beta = 2.0
+        quantity.limiter = "minmod_kurganov"
         quantity.extrapolate_second_order()
 
@@ -405 +407 @@
 
         #Work out the others
-
+        quantity.beta = 2.0
+        quantity.limiter = "minmod_kurganov"
         quantity.extrapolate_second_order()
         
@@ -496 +499 @@
 
 
+
+
+    def test_limit_minmod(self):
+        quantity = Quantity(self.domain4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        #Limit
+        quantity.limiter = 'minmod'
+        quantity.limit_range()
+
+
+
+
+        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
+                                                 [  2.4,   9.6],
+                                                 [ 10.6,  11.4],
+                                                 [ 11.4,  14.6]] )
+
+        from anuga_1d.base.quantity_ext import limit_minmod_ext
+
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        limit_minmod_ext(quantity)
+
+
+
+
+        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
+                                                 [  2.4,   9.6],
+                                                 [ 10.6,  11.4],
+                                                 [ 11.4,  14.6]] )
+
+
+
+    def test_limit_minmod_kurganov(self):
+        quantity = Quantity(self.domain4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        #Limit
+        quantity.limiter = 'minmod_kurganov'
+        quantity.beta = 0.5
+        quantity.limit_range()
+
+
+        #print quantity.vertex_values
+
+        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
+                                                 [  4.2,   7.8],
+                                                 [ 10.8,  11.2],
+                                                 [ 11.4,  14.6]])
+
+        from anuga_1d.base.quantity_ext import limit_minmod_kurganov_ext
+
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        limit_minmod_kurganov_ext(quantity)
+
+
+        #print quantity.vertex_values
+
+
+    def test_limit_vanleer(self):
+        quantity = Quantity(self.domain4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        #Limit
+        quantity.set_limiter('vanleer')
+        quantity.limit_range()
+
+
+        #print quantity.vertex_values
+
+        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
+                                                 [  2.32653061,   9.67346939],
+                                                 [ 10.39393939,  11.60606061],
+                                                 [ 11.4,         14.6       ]] )
+
+        from anuga_1d.base.quantity_ext import limit_vanleer_ext
+
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        limit_vanleer_ext(quantity)
+
+
+        #print quantity.vertex_values
+
+        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
+                                                 [  2.32653061,   9.67346939],
+                                                 [ 10.39393939,  11.60606061],
+                                                 [ 11.4,         14.6       ]] )
+#
 
     def test_find_qmax_qmin(self):
@@ -542 +642 @@
 
         #Extrapolate
+        quantity.beta = 2.0
+        quantity.limiter = "minmod_kurganov"
         quantity.extrapolate_second_order()
 

anuga_work/development/2010-projects/anuga_1d/config.py

@@ -2 +2 @@
 """
 
-epsilon = 1.0e-6
-h0 = 1.0e-6
+epsilon = 1.0e-12
+h0 = 1.0e-12
 
 default_boundary_tag = 'exterior'

anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

@@ -12 +12 @@
 
 h1 = 10.0
-h0 = 0.0
+h0 = 0.1
 
 def analytical_sol(C,t):
@@ -81 +81 @@
 k = 0
 
-N = 800
+N = 3200
 print "Evaluating domain with %d cells" %N
-domain = Domain(*uniform_mesh(N))
+domain = Domain(*uniform_mesh(N,x_1=L))
     
 domain.set_quantity('stage', stage)
@@ -90 +90 @@
 
 domain.set_boundary({'left': Br, 'right' : Br})
-domain.order = 2
-domain.set_timestepping_method('euler')
+domain.set_spatial_order(2)
+domain.set_timestepping_method('rk2')
 domain.set_CFL(1.0)
-domain.set_limiter("minmod")
+domain.set_limiter("minmod_kurganov")
 #domain.h0=0.0001
 
@@ -136 +136 @@
 
 pylab.plot(X.flat,VelocityV.flat)
-plot2.set_ylim([-20,10])
+plot2.set_ylim([-15,15])
 
 pylab.xlabel('Position')

anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -90 +90 @@
         self.quantities_to_be_stored = ['stage','xmomentum']
 
-        self.__doc__ = 'shallow_water_domain'
+        self.__doc__ = 'sww_domain'
 
 
@@ -178 +178 @@
             yield(t)
 
-    def initialise_storage(self):
-        """Create and initialise self.writer object for storing data.
-        Also, save x and bed elevation
-        """
-
-        import data_manager
-
-        #Initialise writer
-        self.writer = data_manager.get_dataobject(self, mode = 'w')
-
-        #Store vertices and connectivity
-        self.writer.store_connectivity()
-
-
-    def store_timestep(self, name):
-        """Store named quantity and time.
-
-        Precondition:
-           self.write has been initialised
-        """
-        self.writer.store_timestep(name)
+
+
 
 
@@ -253 +234 @@
     u_C   = Velocity.centroid_values
         
-    #print id(h_C)
-    #FIXME replace with numpy.where
-    for i in range(N):
-        h_C[i] = w_C[i] - z_C[i]                                                
-        if h_C[i] <= 0.0:
-            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
-            h_C[i] = 1.0e-15
-            w_C[i] = z_C[i]
-            uh_C[i] = 0.0
-                
-    
-##    for i in range(len(h_C)):
-##      if h_C[i] < epsilon:
-##          u_C[i] = 0.0  #Could have been negative
-##          h_C[i] = 0.0
-##      else:
-
-    u_C[:,]  = uh_C/(h_C + h0/h_C)
+    #Calculate height (and fix negatives)better be non-negative!
+    h_C[:] = w_C - z_C
+    u_C[:]  = uh_C/(h_C + h0/h_C)
         
     for name in [ 'velocity', 'stage' ]:
@@ -284 +250 @@
             raise 'Unknown order'
 
+
     stage_V = domain.quantities['stage'].vertex_values                  
     bed_V   = domain.quantities['elevation'].vertex_values      
@@ -291 +258 @@
 
     h_V[:,:]    = stage_V - bed_V
-    for i in range(len(h_C)):
-        for j in range(2):
-            if h_V[i,j] < 0.0 :
-                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
-                dh = h_V[i,j]
-                h_V[i,j] = 0.0
-                stage_V[i,j] = bed_V[i,j]
-                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
-                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
-                
+
+#    # protect from edge values going negative
+#    h_V[:,1] = numpy.where(h_V[:,0] < 0.0 , h_V[:,1]-h_V[:,0], h_V[:,1])
+#    h_V[:,0] = numpy.where(h_V[:,0] < 0.0 , 0.0, h_V[:,0])
+#
+#    h_V[:,0] = numpy.where(h_V[:,1] < 0.0 , h_V[:,0]-h_V[:,1], h_V[:,0])
+#    h_V[:,1] = numpy.where(h_V[:,1] < 0.0 , 0.0, h_V[:,1])
+#
+#
+#    stage_V[:,:] = bed_V + h_V
     xmom_V[:,:] = u_V * h_V
     

anuga_work/development/2010-projects/anuga_1d/utilities/util_ext.h

@@ -31 +31 @@
 
 
+double sign(double x) {
+  //Return sign of a double
+
+  if (x>0.0) return 1.0;
+  else if (x<0.0) return -1.0;
+  else return 0.0;
+}
+
+
 int _gradient(double x0, double y0, 
               double x1, double y1,