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

Timestamp:

May 13, 2013, 5:09:25 PM (12 years ago)

Author:

davies

Message:

Updates to bal_dev

File:

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

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

@@ -66 +66 @@
         self.set_CFL(1.0)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2'#'rk2'#'euler'#'rk2' 
+        self.timestepping_method='rk2'#'rk2'#'rk2'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -74 +74 @@
         self.extrapolate_velocity_second_order=True
 
-        # Note that the extrapolation method used in quantity_ext.c (e.g.
-        # extrapolate_second_order_and_limit_by_edge) uses a constant value for
-        # all the betas.  
         self.beta_w=0.0
         self.beta_w_dry=0.0
@@ -83 +80 @@
         self.beta_vh=0.0
         self.beta_vh_dry=0.0
+
+        #self.epsilon=1.0e-100
 
         #self.optimise_dry_cells=True
@@ -111 +110 @@
         self.boundary_flux_sum=numpy.ndarray(1)
         self.boundary_flux_sum[0]=0.
+
+        self.call=1 # Integer counting how many times we call compute_fluxes_central
+
+        # Integer recording the order of the time-stepping scheme
+        if(self.timestepping_method=='rk2'):
+          self.timestep_order=2
+        elif(self.timestepping_method=='euler'):
+          self.timestep_order=1
+        elif(self.timestepping_method=='rk3'):
+          self.timestep_order=3
+        else:
+          err_mess='ERROR: timestepping_method= ' + self.timestepping_method +' not supported in this solver'
+          raise Exception, err_mess
 
         print '##########################################################################'
@@ -252 +264 @@
         #    raise Exception, err_mess 
 
-        if(domain.timestepping_method=='rk2'):
-          timestep_order=2
-        elif(domain.timestepping_method=='euler'):
-          timestep_order=1
+        #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
+        #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']
@@ -271 +283 @@
         timestep = float(sys.maxint)
 
+        domain.call+=1
         flux_timestep = compute_fluxes_ext(timestep,
                                            domain.epsilon,
@@ -297 +310 @@
                                            domain.max_speed,
                                            int(domain.optimise_dry_cells),
-                                           timestep_order,
+                                           domain.timestep_order,
                                            Stage.centroid_values,
                                            Xmom.centroid_values,
@@ -304 +317 @@
                                            Bed.vertex_values)
 
-        #import pdb
-        #pdb.set_trace()
-        #print 'flux timestep: ', flux_timestep, domain.timestep
+
+        # Update the boundary flux integral
         if(domain.timestepping_method=='rk2'):
-          if(flux_timestep == float(sys.maxint)):
+            if(domain.call%2==1):
               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.
+
+        elif(domain.timestepping_method=='euler'):
+            domain.boundary_flux_integral=0.
+            # This presently doesn't work -- this section of code may need to go elsewhere
+            #domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+            #                                  domain.boundary_flux_sum[0]*domain.timestep
+            #domain.boundary_flux_sum[0]=0.
+
+        elif(domain.timestepping_method=='rk3'):
+            domain.boundary_flux_integral=0.
+            # This needs to be designed
+        else:
+            mess='ERROR: domain.timestepping_method', domain.timestepping_method,' method not recognised'
+            raise Exception, mess
 
         domain.flux_timestep = flux_timestep