Changeset 5175

Timestamp:

Mar 19, 2008, 8:02:05 PM (17 years ago)

Author:

steve

Message:

Ole and I seem to have got edge limiteri working.
Set domain.use_edge_limiter = True
to test the method

Files:

• anuga_core/source/anuga/config.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (5 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (5 diffs)
• anuga_core/source/anuga/shallow_water/test_data_manager.py (modified) (2 diffs)
• anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py (modified) (1 diff)
• anuga_validation/convergence_study/dam_break.py (modified) (3 diffs)
• anuga_validation/convergence_study/wave.py (modified) (2 diffs)
• anuga_work/development/attenuation_near_shore/run_beach_structured.py (modified) (1 diff)

anuga_core/source/anuga/config.py

@@ -176 +176 @@
 optimise_dry_cells = True # Exclude dry and still cells from flux computation
 optimised_gradient_limiter = True # Use hardwired gradient limiter
+use_edge_limiter = False # The edge limiter is better, but most runs have been using vertex limiting 
 
 points_file_block_line_size = 500 # Number of lines read in from a points file

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -104 +104 @@
 from anuga.config import optimise_dry_cells
 from anuga.config import optimised_gradient_limiter
+from anuga.config import use_edge_limiter
 
 #---------------------
@@ -179 +180 @@
 
         # Limiters
-        self.use_old_limiter = True
+        self.use_edge_limiter = use_edge_limiter
 
         self.optimised_gradient_limiter = optimised_gradient_limiter
@@ -395 +396 @@
         # Call correct module function
         # (either from this module or C-extension)
-        if self.use_old_limiter:
-            distribute_to_vertices_and_edges(self)
-        else:
+        if self.use_edge_limiter:            
             for name in self.conserved_quantities:
                 Q = self.quantities[name]
@@ -404 +403 @@
                 elif self._order_ == 2:
                     Q.extrapolate_second_order_and_limit()
-                    if name == 'stage':
-                        Q.bound_vertices_below_by_quantity(self.quantities['elevation'])
                 else:
                     raise 'Unknown order'
- 
+            balance_deep_and_shallow(self)
+
+            #Compute edge values by interpolation
+            for name in self.conserved_quantities:
+                Q = self.quantities[name]
+                Q.interpolate_from_vertices_to_edges()
+        else:
+            distribute_to_vertices_and_edges(self)
+
 
 
@@ -871 +876 @@
     #
     from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
-    
+
+
+    #print 'calling balance depth and shallow'
     # Shortcuts
     wc = domain.quantities['stage'].centroid_values

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -589 +589 @@
     // Calculate minimal depth across all three vertices
     hmin = min(hv[0], min(hv[1], hv[2]));
+
+    //if (hmin < 0.0 ) {
+    //  printf("hmin = %f\n",hmin);
+    //}
     
 
@@ -647 +651 @@
       }
 
-      
+      //printf("alpha=%f, tight_slope_limiters=%d\n", alpha, tight_slope_limiters);
+
       /*      
       // Experimental code for controlling velocities at vertices.
@@ -704 +709 @@
 
         // Update momentum at vertices
-        if (tight_slope_limiters == 1) {        
-          // FIXME(Ole): Here's what I think (as of 17 Nov 2007) 
+
+
+        // Update momentum at vertices
+        // Update momentum as a linear combination of
+        // xmomc and ymomc (shallow) and momentum
+        // from extrapolator xmomv and ymomv (deep).
+
+
+        //xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
+        //ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+
+        if (tight_slope_limiters == 1) {
+          // FIXME(Ole): Here's what I think (as of 17 Nov 2007)
           // we need to do. Simple and efficient:
         
@@ -721 +737 @@
           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
           xmomv[k3+i] = uc*hv[i];
-          ymomv[k3+i] = vc*hv[i];       
+          ymomv[k3+i] = vc*hv[i];
         } else {
           // Update momentum as a linear combination of
           // xmomc and ymomc (shallow) and momentum
           // from extrapolator xmomv and ymomv (deep).
-          // FIXME (Ole): Is this really needed? Could we use the above 
+          // FIXME (Ole): Is this really needed? Could we use the above
           // instead?
           
@@ -1189 +1205 @@
 
       hfactor = 0.0;
-      if (hmin > 0.1 ) {
-          hfactor = (hmin-0.1)/(hmin+0.4);
+      if (hmin > 0.001 ) {
+          hfactor = (hmin-0.001)/(hmin+0.004);
       }
       

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -244 +244 @@
         # Get the variables
         range = fid.variables['stage_range'][:]
+        #print range
         assert allclose(range,[-0.93519, 0.15]) or\
-               allclose(range,[-0.9352743, 0.15]) # Old slope limiters
+               allclose(range,[-0.9352743, 0.15]) or\
+               allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
         
         range = fid.variables['xmomentum_range'][:]
-
-
-
+        #print range
         assert allclose(range,[0,0.4695096]) or \
                allclose(range,[0,0.47790655]) or\
-               allclose(range,[0,0.46941957]) 
+               allclose(range,[0,0.46941957]) or\
+               allclose(range,[0,0.47769409])
+
         
         range = fid.variables['ymomentum_range'][:]
-        #assert allclose(range,[0,0.02174380])
-        
-        assert allclose(range,[0,0.02174439]) or\
+        #print range
+        assert allclose(range,[0,0.02174380]) or\
+               allclose(range,[0,0.02174439]) or\
                allclose(range,[0,0.02283983]) or\
-               allclose(range,[0,0.0217342]) # Old slope limiters 
+               allclose(range,[0,0.0217342]) or\
+               allclose(range,[0,0.0227564]) # Old slope limiters 
         
         fid.close()
@@ -330 +333 @@
 
         extrema = fid.variables['xmomentum.extrema'][:]
-        assert allclose(extrema,[-0.06062178, 0.47873023])
+        assert allclose(extrema,[-0.06062178, 0.47873023]) or allclose(extrema, [-0.06062178, 0.47847986])
         
         extrema = fid.variables['ymomentum.extrema'][:]
-        assert allclose(extrema,[0.00, 0.0625786])        
+        assert allclose(extrema,[0.00, 0.0625786]) or allclose(extrema,[0.00, 0.06062178])
 
         time_interval = fid.variables['extrema.time_interval'][:]

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -2518 +2518 @@
         #print L[1]
         assert allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
-               allclose(L[1], [4.18944138, 16.10955862, 20.001]) # old limiters
+               allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
+               allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
         
         for i in range(len(L)):

anuga_validation/convergence_study/dam_break.py

@@ -41 +41 @@
 # Setup domain
 #------------------------------------------------------------------------------
-dx = 250.
+dx = 200.
 dy = dx
 L = 100000.
-W = dx
+W = 10*dx
 
 # structured mesh
@@ -62 +62 @@
 print domain.get_timestepping_method()
 
-domain.use_old_limiter = True
+domain.use_old_limiter = False
+domain.tight_slope_limiters = True
+
 domain.CFL = 1.0
 
-domain.beta_w      = 1.7
-domain.beta_w_dry  = 0.0
-domain.beta_uh     = 1.7
-domain.beta_uh_dry = 0.0
-domain.beta_vh     = 1.7
-domain.beta_vh_dry = 0.0
-
+domain.beta_w      = 0.6
+domain.beta_uh     = 0.6
+domain.beta_vh     = 0.6
+domain.beta_h    = 0.0
 
 
@@ -120 +119 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 20.0, finaltime = 60*60.):
-    domain.write_time()
+for t in domain.evolve(yieldstep = 100.0, finaltime = 60*60.):
+    print domain.timestepping_statistics(track_speeds=True)
     vis.update()
     

anuga_validation/convergence_study/wave.py

@@ -42 +42 @@
 # Setup domain
 #------------------------------------------------------------------------------
-dx = 500.
+dx = 200.
 dy = dx
 L = 100000.
@@ -58 +58 @@
 # Setup Algorithm
 #------------------------------------------------------------------------------
-#domain.set_timestepping_method('euler')
+domain.set_timestepping_method('rk2')
+domain.set_default_order(2)
 
-#domain.set_default_order(2)
-#domain.set_timestepping_method('euler')
-#domain.set_all_limiters(2.0)
-#domain.use_old_limiter = True
-#domain.CFL = 1.0
+print domain.get_timestepping_method()
 
+domain.use_edge_limiter = True
+domain.tight_slope_limiters = True
 
-domain.set_default_order(2)
-#domain.set_timestepping_method('euler')
+domain.CFL = 1.0
 
-domain.use_old_limiter = True
-
-domain.beta_w      = 2.0
-domain.beta_w_dry  = 0.5
-domain.beta_uh     = 2.0
-domain.beta_uh_dry = 0.5
-domain.beta_vh     = 2.0
-domain.beta_vh_dry = 0.5
+domain.beta_w      = 0.6
+domain.beta_w_dry  = 0.0
+domain.beta_uh     = 0.6
+domain.beta_uh_dry = 0.0
+domain.beta_vh     = 0.6
+domain.beta_vh_dry = 0.0
+domain.beta_h    = 0.0
 
 

anuga_work/development/attenuation_near_shore/run_beach_structured.py

@@ -140 +140 @@
     domain.set_datadir(pro_instance.outputdir)
     domain.set_default_order(2) # Use second order spatial scheme
+
     domain.set_timestepping_method('rk2')
-    domain.beta_w      = 2.0
-    domain.beta_uh     = 2.0
-    domain.beta_vh     = 2.0
-    domain.beta_h      = 2.0
+    domain.set_default_order(2)
+    domain.use_old_limiter = False
+    domain.tight_slope_limiters = True
+    
+    domain.beta_w      = 0.6
+    domain.beta_uh     = 0.6
+    domain.beta_vh     = 0.6
+    domain.beta_h      = 0.0
+
+
     domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     domain.set_minimum_storable_height(0.001)