Changeset 8865 for trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

Timestamp:

May 12, 2013, 8:11:11 PM (12 years ago)

Author:

davies

Message:

Updates to balanced_dev

File:

• trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py (modified) (7 diffs)

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

@@ -59 +59 @@
                             number_of_full_triangles = number_of_full_triangles)
        
+        #self.forcing_terms.append(manning_friction_implicit)
         #------------------------------------------------
         # set some defaults
         # Most of these override the options in config.py
         #------------------------------------------------
-        self.set_CFL(1.00)
+        self.set_CFL(1.0)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2' #'euler'#'rk2'#'euler'#'rk2' 
+        self.timestepping_method='rk2'#'rk2'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -76 +77 @@
         # extrapolate_second_order_and_limit_by_edge) uses a constant value for
         # all the betas.  
-        self.beta_w=1.0
+        self.beta_w=0.0
         self.beta_w_dry=0.0
-        self.beta_uh=1.0
+        self.beta_uh=0.0
         self.beta_uh_dry=0.0
-        self.beta_vh=1.0
+        self.beta_vh=0.0
         self.beta_vh_dry=0.0
 
@@ -241 +242 @@
 
         # Shortcuts
-        if(domain.timestepping_method!='rk2'):
-            err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
-                     'You need to use rk2 timestepping with this solver, ' +\
-                     ' because there is presently a hack in compute fluxes central. \n'+\
-                     ' The HACK: The timestep will only ' +\
-                     'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
-                     ' correct treatment of wetting and drying'
-
-            raise Exception, err_mess 
+        #if(domain.timestepping_method!='rk2'):
+        #    err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
+        #             'You need to use rk2 timestepping with this solver, ' +\
+        #             ' because there is presently a hack in compute fluxes central. \n'+\
+        #             ' The HACK: The timestep will only ' +\
+        #             'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
+        #             ' correct treatment of wetting and drying'
+
+        #    raise Exception, err_mess 
+
+        if(domain.timestepping_method=='rk2'):
+          timestep_order=2
+        elif(domain.timestepping_method=='euler'):
+          timestep_order=1
+        #elif(domain.timestepping_method=='rk3'):
+        #  timestep_order=3
+        else:
+          err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver'
+          raise Exception, err_mess
+
+        #print 'timestep_order=', timestep_order
 
         Stage = domain.quantities['stage']
@@ -260 +273 @@
         flux_timestep = compute_fluxes_ext(timestep,
                                            domain.epsilon,
-                                           domain.H0,
+                                           domain.minimum_allowed_height,
                                            domain.g,
                                            domain.boundary_flux_sum,
@@ -284 +297 @@
                                            domain.max_speed,
                                            int(domain.optimise_dry_cells),
+                                           timestep_order,
                                            Stage.centroid_values,
                                            Xmom.centroid_values,
@@ -293 +307 @@
         #pdb.set_trace()
         #print 'flux timestep: ', flux_timestep, domain.timestep
-        if(flux_timestep == float(sys.maxint)):
-            domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
-                                              domain.boundary_flux_sum[0]*domain.timestep*0.5
-            #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
-            domain.boundary_flux_sum[0]=0.
+        if(domain.timestepping_method=='rk2'):
+          if(flux_timestep == float(sys.maxint)):
+              domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+                                                domain.boundary_flux_sum[0]*domain.timestep*0.5
+              #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
+              domain.boundary_flux_sum[0]=0.
 
         domain.flux_timestep = flux_timestep
@@ -590 +605 @@
 
 
+#def manning_friction_implicit(domain):
+#    """Apply (Manning) friction to water momentum
+#    Wrapper for c version
+#    """
+#
+#    from shallow_water_ext import manning_friction_flat
+#    from shallow_water_ext import manning_friction_sloped
+#
+#    xmom = domain.quantities['xmomentum']
+#    ymom = domain.quantities['ymomentum']
+#
+#    x = domain.get_vertex_coordinates()
+#    
+#    w = domain.quantities['stage'].centroid_values
+#    z = domain.quantities['elevation'].vertex_values
+#
+#    uh = xmom.centroid_values
+#    vh = ymom.centroid_values
+#    eta = domain.quantities['friction'].centroid_values
+#
+#    xmom_update = xmom.semi_implicit_update
+#    ymom_update = ymom.semi_implicit_update
+#
+#    eps = domain.epsilon
+#    g = domain.g
+#
+#    if domain.use_sloped_mannings:
+#        manning_friction_sloped(g, eps, x, w, uh, vh, z, eta, xmom_update, \
+#                                ymom_update)
+#    else:
+#        manning_friction_flat(g, eps, w, uh, vh, z, eta, xmom_update, \
+#                                #ymom_update)