Changeset 4677 for anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

Timestamp:

Aug 22, 2007, 4:39:07 PM (17 years ago)

Author:

ole

Message:

Installed improvement to flux limiter ensuring consistency between
velocity and momentum (not ensured before).

Also implemented optional policy for treating isolated instances of
very small timesteps (this is experimental and disabled by default).

File:

• anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py (modified) (5 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -171 +171 @@
 
         #Defaults
-        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon, CFL
+        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
+        from anuga.config import CFL
+        from anuga.config import protect_against_isolated_degenerate_timesteps
         self.beta_w = beta_w
         self.beta_h = beta_h
         self.epsilon = epsilon
+        self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
+        
 
         #FIXME: Maybe have separate orders for h-limiter and w-limiter?
@@ -646 +650 @@
             for name in self.quantities:
                 q = self.quantities[name]
-                X,Y,A,V = q.get_vertex_values()                
+
+                V = q.get_values(location='vertices', indices=[k])[0]
+                C = q.get_values(location='centroids', indices=[k])                 
                 
                 s = '    %s:\t %.4f,\t %.4f,\t %.4f,\t %.4f\n'\
-                    %(name, A[3*k], A[3*k+1], A[3*k+2], q.get_values(location='centroids')[k]) 
+                    %(name, V[0], V[1], V[2], C[0])                 
 
                 msg += s
@@ -962 +968 @@
         from anuga.config import min_timestep, max_timestep
 
+        
+        
+        # Protect against degenerate timesteps arising from isolated
+        # triangles
+        if self.protect_against_isolated_degenerate_timesteps is True and\
+               self.max_speed > 10.0:
+
+            # Setup 10 bins for speed histogram
+            from anuga.utilities.numerical_tools import histogram, create_bins
+        
+            bins = create_bins(self.max_speed, 10)
+            hist = histogram(self.max_speed, bins)
+
+            # Look for characteristic signature
+            if len(hist) > 1 and\
+                hist[-1] > 0 and\
+                hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
+                    # Danger of isolated degenerate triangles
+                    # print self.timestepping_statistics(track_speeds=True)  
+                    
+                    # Find triangles in last bin
+                    # FIXME - speed up using Numeric
+                    d = 0
+                    for i in range(self.number_of_full_triangles):
+                        if self.max_speed[i] > bins[-1]:
+                            msg = 'Time=%f: Ignoring isolated high speed triangle ' %self.time
+                            msg += '#%d of %d with max speed=%f'\
+                                  %(i, self.number_of_full_triangles, self.max_speed[i])
+                            #print msg
+                    
+                            # print 'Found offending triangle', i, self.max_speed[i] 
+                            self.get_quantity('xmomentum').set_values(0.0, indices=[i])
+                            self.get_quantity('ymomentum').set_values(0.0, indices=[i])
+                            self.max_speed[i]=0.0
+                            d += 1
+                    
+                    #print 'Adjusted %d triangles' %d        
+                    #print self.timestepping_statistics(track_speeds=True)      
+                
+
+                        
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
@@ -970 +1017 @@
         self.min_timestep = min(timestep, self.min_timestep)
 
+           
+ 
         #Protect against degenerate time steps
         if timestep < min_timestep:
@@ -988 +1037 @@
 
 
-                    print self.timestepping_statistics(track_speeds=True)
-
-
-                    raise Exception
+                    #print self.timestepping_statistics(track_speeds=True)
+
+
+                    raise Exception, msg
                 else:
                     #Try to overcome situation by switching to 1 order