Changeset 8404

Timestamp:

Apr 18, 2012, 4:51:48 PM (13 years ago)

Author:

steve

Message:

added in a variable to set higher/lower order algorithm. Also set manning sloped as default

Location:

trunk/anuga_core/source/anuga

Files:

• config.py (modified) (4 diffs)
• file_conversion/asc2dem.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (6 diffs)
• shallow_water/test_forcing.py (modified) (1 diff)
• shallow_water/test_shallow_water_domain.py (modified) (4 diffs)
• shallow_water/test_sww_interrogate.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/config.py

@@ -55 +55 @@
 # values from centroid values
 ################################################################################
-
-# Betas [0;1] control the allowed steepness of gradient for second order
-# extrapolations. Values of 1 allow the steepes gradients while
-# lower values are more conservative. Values of 0 correspond to
-# 1'st order extrapolations.
-#
-
-# There are separate betas for the w, uh, and vh limiters
-# I think these are better SR but they conflict with the unit tests!
-beta_w      = 1.0
-beta_w_dry  = 0.2
-beta_uh     = 1.0
-beta_uh_dry = 0.2
-beta_vh     = 1.0
-beta_vh_dry = 0.2
+# Note the individual beta values are set in domain.set_flow_method which also sets
+# the timestepping method
+
+beta_w = 1.0
 
 # Alpha_balance controls how limiters are balanced between deep and shallow.
@@ -113 +102 @@
 ################################################################################
 
-sloped_mannings_function = False
+sloped_mannings_function = True
 
 ################################################################################
@@ -121 +110 @@
 CFL = 1.0  # CFL condition assigned to domain.CFL - controls timestep size
       
-# Choose type of timestepping,
-#timestepping_method = 'rk2'   # 2nd Order TVD scheme
-timestepping_method = 'euler' # 1st order euler
+# Choose type of timestepping and spatial reconstruction method
+
+timestepping_method = 1
+
+# For shallow water we have a method that sets both timestepping and spatial reconstruction and
+# beta values. In this case the settings for timestepping_method will be overriden 
+
+#flow_algorithm = 1   # 1st order euler and conservative piecewise constant spatial reconstruction
+flow_algorithm = 1.5 # 1st order euler and conservative piecewise linear spatial reconstruction
+#flow_algorithm = 2   # 2nd order TVD scheme and more aggressive piecewise linear spatial reconstruction
+#flow_algorithm = 2.5 # 3rd order TVD scheme and more aggressive piecewise linear spatial reconstruction
+
+
+
 
 # rk2 is a little more stable than euler, so rk2 timestepping
@@ -191 +191 @@
 ################################################################################
 
-use_psyco = True      # Use psyco optimisations
+use_psyco = False      # Use psyco optimisations
 
 optimise_dry_cells = True # Exclude dry and still cells from flux computation

trunk/anuga_core/source/anuga/file_conversion/asc2dem.py

@@ -147 +147 @@
         raise Exception, msg
 
-    NODATA_value = int(lines[5].split()[1].strip())
+    NODATA_value = int(float(lines[5].split()[1].strip()))
 
     assert len(lines) == nrows + 6

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -217 +217 @@
         from anuga.config import minimum_storable_height
         from anuga.config import minimum_allowed_height, maximum_allowed_speed
-        from anuga.config import g, beta_w, beta_w_dry, \
-             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
+        from anuga.config import g
+        from anuga.config import tight_slope_limiters
         from anuga.config import extrapolate_velocity_second_order
         from anuga.config import alpha_balance
@@ -227 +227 @@
         from anuga.config import compute_fluxes_method
         from anuga.config import sloped_mannings_function
+        from anuga.config import flow_algorithm
         
 
@@ -232 +233 @@
         self.maximum_allowed_speed = maximum_allowed_speed
 
-        self.set_extrapolate_velocity(extrapolate_velocity_second_order)
-
+
+
+        self.minimum_storable_height = minimum_storable_height
         self.g = g
-        self.beta_w = beta_w
-        self.beta_w_dry = beta_w_dry
-        self.beta_uh = beta_uh
-        self.beta_uh_dry = beta_uh_dry
-        self.beta_vh = beta_vh
-        self.beta_vh_dry = beta_vh_dry
+
         self.alpha_balance = alpha_balance
-
         self.tight_slope_limiters = tight_slope_limiters
         self.optimise_dry_cells = int(optimise_dry_cells)
-
-
-        self.minimum_storable_height = minimum_storable_height
-
-         # Limiters
+        self.set_extrapolate_velocity(extrapolate_velocity_second_order)
+
+
         self.set_use_edge_limiter(use_edge_limiter)
         self.optimised_gradient_limiter = optimised_gradient_limiter
@@ -255 +249 @@
 
         self.set_sloped_mannings_function(sloped_mannings_function)
-
         self.set_compute_fluxes_method(compute_fluxes_method)
+        self.set_flow_algorithm(flow_algorithm)
+
+
+    def get_algorithm_parameters(self):
+        """Get the standard parameter that arecurently set (as a dictionary)
+        """
+
+        parameters = {}
+
+        parameters['minimum_allowed_height']  = self.minimum_allowed_height
+        parameters['maximum_allowed_speed']   = self.maximum_allowed_speed
+        parameters['minimum_storable_height'] = self.minimum_storable_height
+        parameters['g']                       = self.g
+        parameters['alpha_balance']           = self.alpha_balance
+        parameters['tight_slope_limiters']    = self.tight_slope_limiters
+        parameters['optimise_dry_cells']      = self.optimise_dry_cells
+        parameters['use_edge_limiter']        = self.use_edge_limiter
+        parameters['use_centroid_velocities'] = self.use_centroid_velocities
+        parameters['use_sloped_mannings']     = self.use_sloped_mannings
+        parameters['compute_fluxes_method']   = self.get_compute_fluxes_method()
+        parameters['flow_algorithm']             = self.get_flow_algorithm()
+
+        parameters['optimised_gradient_limiter']        = self.optimised_gradient_limiter
+        parameters['extrapolate_velocity_second_order'] = self.extrapolate_velocity_second_order
+        
+        return parameters
+
+    def print_algorithm_parameters(self):
+        """Print the standard parameters that are curently set (as a dictionary)
+        """
+
+        print '#============================'
+        print '# Domain Algorithm Parameters '
+        print '#============================'
+        from pprint import pprint
+        pprint(self.get_algorithm_parameters(),indent=4)
+
+        print '#----------------------------'
 
 
@@ -310 +341 @@
 
         Currently
-           'original'
-           'well_balanced_1
+           original
+           wb_1
         """
 
         return self.compute_fluxes_method
 
+
+    def set_flow_algorithm(self, flag=1.5):
+        """Set combination of slope limiting and time stepping
+
+        Currently
+           1
+           1.5
+           2
+           2.5
+        """
+
+        flag = str(flag)
+
+        flow_algorithms = ['1', '1.5', '2', '2.5']
+
+        if flag in flow_algorithms:
+            self.flow_algorithm = flag
+        else:
+            msg = 'Unknown flow_algorithm. \nPossible choices are:\n'+ \
+            ', '.join(flow_algorithm)+'.'
+            raise Exception(msg)
+
+
+        if self.flow_algorithm == '1':
+            self.set_timestepping_method(1)
+            self.set_default_order(1)
+
+        if self.flow_algorithm == '1.5':
+            self.set_timestepping_method(1)
+            self.set_default_order(2)
+            beta_w      = 1.0
+            beta_w_dry  = 0.2
+            beta_uh     = 1.0
+            beta_uh_dry = 0.2
+            beta_vh     = 1.0
+            beta_vh_dry = 0.2
+            self.set_betas(beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry)
+
+
+        if self.flow_algorithm == '2':
+            self.set_timestepping_method(2)
+            self.set_default_order(2)
+            beta_w      = 1.5
+            beta_w_dry  = 0.2
+            beta_uh     = 1.5
+            beta_uh_dry = 0.2
+            beta_vh     = 1.5
+            beta_vh_dry = 0.2
+            self.set_betas(beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry)
+
+        if self.flow_algorithm == '2.5':
+            self.set_timestepping_method(3)
+            self.set_default_order(2)
+            beta_w      = 1.5
+            beta_w_dry  = 0.2
+            beta_uh     = 1.5
+            beta_uh_dry = 0.2
+            beta_vh     = 1.5
+            beta_vh_dry = 0.2
+            self.set_betas(beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry)
+
+
+
+    def get_flow_algorithm(self):
+        """Get method used for timestepping and spatial discretisation
+
+        Currently
+           1, 1.5, 2, 2.5
+        """
+
+        return self.flow_algorithm
+    
 
     def set_gravity_method(self):
@@ -393 +496 @@
         self.beta_vh_dry = beta
         self.quantities['ymomentum'].beta = beta
+
+
+    def set_betas(self, beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry):
+        """Assign beta values in the range  [0,2] to all limiters.
+        0 Corresponds to first order, where as larger values make use of
+        the second order scheme.
+        """
+
+        self.beta_w = beta_w
+        self.beta_w_dry = beta_w_dry
+        self.quantities['stage'].beta = beta_w
+
+        self.beta_uh = beta_uh
+        self.beta_uh_dry = beta_uh_dry
+        self.quantities['xmomentum'].beta = beta_uh
+
+        self.beta_vh = beta_vh
+        self.beta_vh_dry = beta_vh_dry
+        self.quantities['ymomentum'].beta = beta_vh
+
+
 
 

trunk/anuga_core/source/anuga/shallow_water/test_forcing.py

@@ -1982 +1982 @@
         domain = Domain(points, vertices)
 
+        # Use the old function which doesn't take into account the extra
+        # wetted area due to slope of bed
+        domain.set_sloped_mannings_function(False)
+
         B = Reflective_boundary(domain)
         domain.set_boundary( {'exterior': B})

trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -409 +409 @@
         assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
 
+
+
+    def test_algorithm_parameters(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
+        #msg = 'The class %s is not a subclass of the generic domain class %s'\
+        #      %(DomainClass, Domain)
+        #assert issubclass(DomainClass, Domain), msg
+
+        domain = Domain(points, vertices)
+        domain.check_integrity()
+
+        for name in ['stage', 'xmomentum', 'ymomentum',
+                     'elevation', 'friction']:
+            assert domain.quantities.has_key(name)
+
+
+
+        parameters = domain.get_algorithm_parameters()
+
+
+
+
+        assert parameters['minimum_allowed_height']  == domain.minimum_allowed_height
+        assert parameters['maximum_allowed_speed']   == domain.maximum_allowed_speed
+        assert parameters['minimum_storable_height'] == domain.minimum_storable_height
+        assert parameters['g']                       == domain.g
+        assert parameters['alpha_balance']           == domain.alpha_balance
+        assert parameters['tight_slope_limiters']    == domain.tight_slope_limiters
+        assert parameters['optimise_dry_cells']      == domain.optimise_dry_cells
+        assert parameters['use_edge_limiter']        == domain.use_edge_limiter
+        assert parameters['use_centroid_velocities'] == domain.use_centroid_velocities
+        assert parameters['use_sloped_mannings']     == domain.use_sloped_mannings
+        assert parameters['compute_fluxes_method']   == domain.get_compute_fluxes_method()
+        assert parameters['flow_algorithm']          == domain.get_flow_algorithm()
+
+        assert parameters['optimised_gradient_limiter']        == domain.optimised_gradient_limiter
+        assert parameters['extrapolate_velocity_second_order'] == domain.extrapolate_velocity_second_order
+
+
+
+
     def test_boundary_conditions(self):
         a = [0.0, 0.0]
@@ -2670 +2722 @@
 
 
-        print domain.quantities['xmomentum'].explicit_update
-        print domain.quantities['ymomentum'].explicit_update
+        #print domain.quantities['xmomentum'].explicit_update
+        #print domain.quantities['ymomentum'].explicit_update
 
         assert num.allclose(domain.quantities['stage'].explicit_update, 0)
         assert num.allclose(domain.quantities['xmomentum'].explicit_update,
-                                [ -2.94,  -2.94, -10.29, -10.29])
+                                [ -2.94,  -2.94, -2.94, -2.94])
         assert num.allclose(domain.quantities['ymomentum'].explicit_update,
-                                [  7.35, 0.0, 0.0, -7.35])
+                                [  0.0, 0.0, 0.0, 0.0])
 
 
@@ -5232 +5284 @@
             pass
 
+#        print domain.quantities['stage'].centroid_values[:4]
+#        print domain.quantities['xmomentum'].centroid_values[:4]
+#        print domain.quantities['ymomentum'].centroid_values[:4]
         assert num.allclose(domain.quantities['stage'].centroid_values[:4],
-                            [0.00206836, 0.01296714, 0.00363415, 0.01438924])
+                            [ 0.001362,    0.01344294,  0.00308829, 0.01470289])
         assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
-                            [0.01360154, 0.00671133, 0.01264578, 0.00648503])
+                            [ 0.01300239,  0.00537933,  0.01214676,  0.00515825])
         assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
-                            [-1.19201077e-003, -7.23647546e-004,
-                             -6.39083123e-005, 6.29815168e-005])
+                            [ -1.13165691e-03,  -6.55330189e-04,
+                              -6.62804076e-05,   5.26313051e-05])
+
+        # old values pre revision 8402
+#        assert num.allclose(domain.quantities['stage'].centroid_values[:4],
+#                            [0.00206836, 0.01296714, 0.00363415, 0.01438924])
+#        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
+#                            [0.01360154, 0.00671133, 0.01264578, 0.00648503])
+#        assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
+#                            [-1.19201077e-003, -7.23647546e-004,
+#                             -6.39083123e-005, 6.29815168e-005])
 
         os.remove(domain.get_name() + '.sww')
@@ -7904 +7968 @@
 if __name__ == "__main__":
     #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
-    suite = unittest.makeSuite(Test_Shallow_Water, 'test_compute_fluxes_')
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test_')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

trunk/anuga_core/source/anuga/shallow_water/test_sww_interrogate.py

@@ -85 +85 @@
         location = get_maximum_inundation_location(swwfile)
         #print 'Runup, location', runup, location
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15) or num.allclose(location[0], 10)
+        assert num.allclose(runup, 6) or num.allclose(runup, 12) # old limiters
+        assert num.allclose(location[0], 40.0) or num.allclose(location[0], 10)
 
         # Check final runup
@@ -153 +153 @@
         runup = get_maximum_inundation_elevation(swwfile)
         location = get_maximum_inundation_location(swwfile)
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15+E) or num.allclose(location[0], 10+E)
+
+        #print runup, location
+        
+        assert num.allclose(runup, 6) or num.allclose(runup, 12) # old limiters
+        assert num.allclose(location[0], 40+E) or num.allclose(location[0], 10+E)
 
         # Check final runup