Changeset 8381

Timestamp:

Apr 3, 2012, 7:45:42 AM (13 years ago)

Author:

steve

Message:

Reorganised the compute_fluxes to incorporate teh gravity terms, and to
have the ability to change method. 3 methods now. ome of the unit test are failing
but will fix those up once I get to work. The code seems to be working ok

Location:

trunk

Files:

• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (10 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (2 diffs)
• anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py (modified) (6 diffs)
• anuga_work/shallow_water_balanced_steve/run_sw_merimbula.py (modified) (1 diff)

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

@@ -140 +140 @@
             other_quantities = ['elevation', 'friction', 'height',
                                 'xvelocity', 'yvelocity', 'x', 'y']
+
+
         
         
@@ -167 +169 @@
 
         #-------------------------------
-        # Operators and Forcing Terms
+        # Forcing Terms
+        #
+        # Gravity is now incorporated in
+        # compute_fluxes routine
         #-------------------------------
         self.forcing_terms.append(manning_friction_implicit)
-        #self.forcing_terms.append(gravity_wb)
-        self.forcing_terms.append(gravity)
-
-
-
+
+
+        #-------------------------------
+        # Operators
+        #-------------------------------
         self.fractional_step_operators = []
         self.kv_operator = None
@@ -250 +255 @@
         self.use_centroid_velocities = use_centroid_velocities
 
+
         self.set_compute_fluxes_method(compute_fluxes_method)
 
@@ -285 +291 @@
            'well_balanced_1
         """
-        compute_fluxes_methods = ['original', 'well_balanced_1']
+        compute_fluxes_methods = ['original', 'wb_1', 'wb_2']
 
         if flag in compute_fluxes_methods:
@@ -295 +301 @@
 
 
+
+
+
     def get_compute_fluxes_method(self):
         """Get method for computing fluxes.
@@ -304 +313 @@
 
         return self.compute_fluxes_method
+
+
+    def set_gravity_method(self):
+        """Gravity method is determined by the compute_fluxes_method
+        """
+
+        if  self.get_compute_fluxes_method() == 'original':
+            self.forcing_terms[0] = gravity
+
+        elif self.get_compute_fluxes_method() == 'wb_1':
+            self.forcing_terms[0] = gravity_wb
+
+        elif self.get_compute_fluxes_method() == 'wb_2':
+            self.forcing_terms[0] = gravity
+
+        else:
+            raise Exception('undefined compute_fluxes method')
 
 
@@ -625 +651 @@
         if self.compute_fluxes_method == 'original':
             from shallow_water_ext import compute_fluxes_ext_central_structure
+            from shallow_water_ext import gravity as gravity_c
+
             self.flux_timestep = compute_fluxes_ext_central_structure(self)
-        elif self.compute_fluxes_method == 'well_balanced_1':
+            gravity_c(self)
+
+
+        elif self.compute_fluxes_method == 'wb_1':
             from shallow_water_ext import compute_fluxes_ext_wb
+            from shallow_water_ext import gravity as gravity_c
+
             self.flux_timestep = compute_fluxes_ext_wb(self)
+            gravity_c(self)
+
+        elif self.compute_fluxes_method == 'wb_2':
+            from shallow_water_ext import compute_fluxes_ext_central_structure
+            from shallow_water_ext import gravity_wb as gravity_wb_c
+
+            self.flux_timestep = compute_fluxes_ext_central_structure(self)
+            gravity_wb_c(self)
+
         else:
             raise Exception('unknown compute_fluxes_method')
+
+
 
 
@@ -1087 +1131 @@
         
         return message        
-           
+          
+
+
 ################################################################################
 # End of class Shallow Water Domain
@@ -1096 +1142 @@
 #-----------------
 
-def compute_fluxes(domain):
-    """Compute fluxes and timestep suitable for all volumes in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
-
-    This wrapper calls the underlying C version of compute fluxes
-    """
-
-    import sys
-    from shallow_water_ext import compute_fluxes_ext_central \
-                                  as compute_fluxes_ext
-
-    # Shortcuts
-    Stage = domain.quantities['stage']
-    Xmom = domain.quantities['xmomentum']
-    Ymom = domain.quantities['ymomentum']
-    Bed = domain.quantities['elevation']
-
-
-
-    timestep = float(sys.maxint)
-
-    flux_timestep = compute_fluxes_ext(timestep,
-                                       domain.epsilon,
-                                       domain.H0,
-                                       domain.g,
-                                       domain.neighbours,
-                                       domain.neighbour_edges,
-                                       domain.normals,
-                                       domain.edgelengths,
-                                       domain.radii,
-                                       domain.areas,
-                                       domain.tri_full_flag,
-                                       Stage.edge_values,
-                                       Xmom.edge_values,
-                                       Ymom.edge_values,
-                                       Bed.edge_values,
-                                       Stage.boundary_values,
-                                       Xmom.boundary_values,
-                                       Ymom.boundary_values,
-                                       Stage.explicit_update,
-                                       Xmom.explicit_update,
-                                       Ymom.explicit_update,
-                                       domain.already_computed_flux,
-                                       domain.max_speed,
-                                       domain.optimise_dry_cells)
-
-    domain.flux_timestep = flux_timestep
-
-
-
-def compute_fluxes_structure(domain):
-    """Compute fluxes and timestep suitable for all volumes in domain.
-
-    Compute total flux for each conserved quantity using "flux_function"
-
-    Fluxes across each edge are scaled by edgelengths and summed up
-    Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities
-    stage, xmomentum and ymomentum
-
-    The maximal allowable speed computed by the flux_function for each volume
-    is converted to a timestep that must not be exceeded. The minimum of
-    those is computed as the next overall timestep.
-
-    Post conditions:
-      domain.explicit_update is reset to computed flux values
-      domain.flux_timestep is set to the largest step satisfying all volumes.
-
-    This wrapper calls the underlying C version of compute fluxes
-    """
-
-
-    from shallow_water_ext import compute_fluxes_ext_central_structure
-
-
-    domain.flux_timestep = compute_fluxes_ext_central_structure(domain)
-
+#def compute_fluxes(domain):
+#    """Compute fluxes and timestep suitable for all volumes in domain.
+#
+#    Compute total flux for each conserved quantity using "flux_function"
+#
+#    Fluxes across each edge are scaled by edgelengths and summed up
+#    Resulting flux is then scaled by area and stored in
+#    explicit_update for each of the three conserved quantities
+#    stage, xmomentum and ymomentum
+#
+#    The maximal allowable speed computed by the flux_function for each volume
+#    is converted to a timestep that must not be exceeded. The minimum of
+#    those is computed as the next overall timestep.
+#
+#    Post conditions:
+#      domain.explicit_update is reset to computed flux values
+#      domain.timestep is set to the largest step satisfying all volumes.
+#
+#    This wrapper calls the underlying C version of compute fluxes
+#    """
+#
+#    import sys
+#    from shallow_water_ext import compute_fluxes_ext_central \
+#                                  as compute_fluxes_ext
+#
+#    # Shortcuts
+#    Stage = domain.quantities['stage']
+#    Xmom = domain.quantities['xmomentum']
+#    Ymom = domain.quantities['ymomentum']
+#    Bed = domain.quantities['elevation']
+#
+#
+#
+#    timestep = float(sys.maxint)
+#
+#    flux_timestep = compute_fluxes_ext(timestep,
+#                                       domain.epsilon,
+#                                       domain.H0,
+#                                       domain.g,
+#                                       domain.neighbours,
+#                                       domain.neighbour_edges,
+#                                       domain.normals,
+#                                       domain.edgelengths,
+#                                       domain.radii,
+#                                       domain.areas,
+#                                       domain.tri_full_flag,
+#                                       Stage.edge_values,
+#                                       Xmom.edge_values,
+#                                       Ymom.edge_values,
+#                                       Bed.edge_values,
+#                                       Stage.boundary_values,
+#                                       Xmom.boundary_values,
+#                                       Ymom.boundary_values,
+#                                       Stage.explicit_update,
+#                                       Xmom.explicit_update,
+#                                       Ymom.explicit_update,
+#                                       domain.already_computed_flux,
+#                                       domain.max_speed,
+#                                       domain.optimise_dry_cells)
+#
+#    domain.flux_timestep = flux_timestep
+#
+#
+#
+#def compute_fluxes_structure(domain):
+#    """Compute fluxes and timestep suitable for all volumes in domain.
+#
+#    Compute total flux for each conserved quantity using "flux_function"
+#
+#    Fluxes across each edge are scaled by edgelengths and summed up
+#    Resulting flux is then scaled by area and stored in
+#    explicit_update for each of the three conserved quantities
+#    stage, xmomentum and ymomentum
+#
+#    The maximal allowable speed computed by the flux_function for each volume
+#    is converted to a timestep that must not be exceeded. The minimum of
+#    those is computed as the next overall timestep.
+#
+#    Post conditions:
+#      domain.explicit_update is reset to computed flux values
+#      domain.flux_timestep is set to the largest step satisfying all volumes.
+#
+#    This wrapper calls the underlying C version of compute fluxes
+#    """
+#
+#
+#    from shallow_water_ext import compute_fluxes_ext_central_structure
+#
+#
+#    domain.flux_timestep = compute_fluxes_ext_central_structure(domain)
+#
 
 
@@ -1401 +1447 @@
 
 
-def gravity_wb(domain):
-    """Apply gravitational pull in the presence of bed slope
-    Wrapper calls underlying C implementation
-    """
-
-    from shallow_water_ext import gravity_wb as gravity_wb_c
-
-    gravity_wb_c(domain)
-
-
-def gravity(domain):
-    """Apply gravitational pull in the presence of bed slope
-    Wrapper calls underlying C implementation
-    """
-
-    from shallow_water_ext import gravity as gravity_c
-
-    gravity_c(domain)
 
 def manning_friction_implicit(domain):

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1601 +1601 @@
 
 
-        printf("h0,1,2 %f %f %f\n",hh[0],hh[1],hh[2]);
+        //printf("h0,1,2 %f %f %f\n",hh[0],hh[1],hh[2]);
         
         // Calculate the side correction term
@@ -1610 +1610 @@
             n1 = D.normals[k6+2*i+1];
 
-            printf("n0, n1 %i %g %g\n",i,n0,n1);
-            fact = 0.5*g*hh[i]*hh[i]*D.edgelengths[k3+i];
+            //printf("n0, n1 %i %g %g\n",i,n0,n1);
+            fact = -0.5*g*hh[i]*hh[i]*D.edgelengths[k3+i];
             sidex = sidex + fact*n0;
             sidey = sidey + fact*n1;

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

@@ -2167 +2167 @@
             assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
 
-        domain.compute_forcing_terms()
-        #from shallow_water_ext import gravity
-        #gravity(domain)
+        #domain.compute_forcing_terms()
+        from shallow_water_ext import gravity
+        gravity(domain)
 
         #print domain.quantities['xmomentum'].explicit_update
@@ -2181 +2181 @@
 
 
-
-    def Xtest_gravity_wb(self):
+    def test_gravity_2(self):
         #Assuming no friction
 
@@ -2204 +2203 @@
             return 3*x
 
-        h = 0.1
+        h = 15
         def stage(x, y):
-            return slope(x,y)+h
+            return h
 
         domain.set_quantity('elevation', slope)
@@ -2215 +2214 @@
             assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
 
-        from shallow_water_ext import gravity_wb
-        gravity_wb(domain)
-
-
-        print domain.quantities['xmomentum'].explicit_update
-        print domain.quantities['ymomentum'].explicit_update
+        #domain.compute_forcing_terms()
+        from shallow_water_ext import gravity
+        gravity(domain)
+
+        #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,
-                            -g*h*3)
+                            [-382.2, -323.4, -205.8, -382.2])
         assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
 
-
-    def Xtest_gravity_wb_2(self):
-        #Assuming no friction
+    def test_gravity_3(self):
+        #Showing standard gravity term is not well balanced
 
         from anuga.config import g
@@ -2249 +2249 @@
             return 3*x
 
+        h = 15
+        def stage(x, y):
+            return h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+        for name in domain.conserved_quantities:
+            assert num.allclose(domain.quantities[name].explicit_update, 0)
+            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
+
+        domain.set_compute_fluxes_method('original')
+        domain.compute_forcing_terms()
+
+
+        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,
+                            [-382.2, -323.4, -205.8, -382.2])
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
+
+    def test_gravity_wb(self):
+        #Assuming no friction
+
+        from anuga.config import g
+
+        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]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
+        h = 0.1
+        def stage(x, y):
+            return slope(x,y)+h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+        for name in domain.conserved_quantities:
+            assert num.allclose(domain.quantities[name].explicit_update, 0)
+            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
+
+        from shallow_water_ext import gravity_wb
+        gravity_wb(domain)
+
+
+        #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,
+                            -g*h*3)
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
+
+    def test_gravity_wb_2(self):
+        #Assuming no friction
+
+        from anuga.config import g
+
+        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]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
         h = 15.0
         def stage(x, y):
@@ -2264 +2358 @@
 
 
-        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,
-                            -g*h*3)
+                            [-396.9, -308.7, -220.5, -396.9])
         assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
 

trunk/anuga_work/shallow_water_balanced_steve/run_sw_merimbula.py

@@ -97 +97 @@
 #domain.set_CFL(0.7)
 #domain.set_beta(1.5)
-domain.set_compute_fluxes_method('original')
+domain.set_compute_fluxes_method('wb_2')
 domain.set_name('merimbula')