Changeset 4713

Timestamp:

Sep 10, 2007, 4:09:23 PM (17 years ago)

Author:

steve

Message:

Implemented 2nd and 3rd order timestepping via the domain.set_timestepping_method method

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py (modified) (9 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (2 diffs)
• anuga_work/development/demos/netherlands.py (modified) (4 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -208 +208 @@
         self.min_timestep = self.max_timestep = 0.0
         self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        self.timestep = 0.0
+        self.flux_timestep = 0.0
 
         # Monitoring
@@ -242 +244 @@
 
 
-    def set_default_order(self, n):
-        """Set default (spatial) order to either 1 or 2
-        """
-
-        msg = 'Default order must be either 1 or 2. I got %s' %n
-        assert n in [1,2], msg
-
-        self.default_order = n
-        self._order_ = self.default_order
+
 
 
@@ -295 +289 @@
 
         return self.time
+
+    def set_default_order(self, n):
+        """Set default (spatial) order to either 1 or 2
+        """
+
+        msg = 'Default order must be either 1 or 2. I got %s' %n
+        assert n in [1,2], msg
+
+        self.default_order = n
+        self._order_ = self.default_order
+        
 
     def set_quantity_vertices_dict(self, quantity_dict):
@@ -964 +969 @@
     def set_timestepping_method(self,timestepping_method):
         
-        if timestepping_method in ['euler', 'rk2']:
+        if timestepping_method in ['euler', 'rk2', 'rk3']:
             self.timestepping_method = timestepping_method
             return
@@ -1113 +1118 @@
                 
             elif self.get_timestepping_method() == 'rk2':
-                self.evolve_one_rk2_step(yieldstep,finaltime)               
-
+                self.evolve_one_rk2_step(yieldstep,finaltime)
+
+            elif self.get_timestepping_method() == 'rk3':
+                self.evolve_one_rk3_step(yieldstep,finaltime)               
+            
             # Update extrema if necessary (for reporting)
             self.update_extrema()            
@@ -1183 +1191 @@
 
 
-    def evolve_one_rk2_step(self,yieldstep, finaltime):
+    def evolve_one_rk2_step(self, yieldstep, finaltime):
         """One 2nd order RK timestep"""
 
@@ -1226 +1234 @@
         #------------------------------------
         #Combine initial and final values
-        #of conserved quantities
+        #of conserved quantities and cleanup
         #------------------------------------
-
+        #combine steps
         self.saxpy_conserved_quantities(0.5, 0.5)
+ 
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+
+
+    def evolve_one_rk3_step(self, yieldstep, finaltime):
+        """One 2nd order RK timestep"""
+
+
+
+        #Save initial initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        initial_time = self.time
+        
+        #--------------------------------------
+        #First euler step
+        #--------------------------------------
+
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+        #Update time
+        self.time += self.timestep
 
         #------------------------------------
-        #Clean up rk step
+        #Second Euler step
         #------------------------------------
             
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine final and initial values
+        #of conserved quantities and cleanup
+        #------------------------------------
+        #combine steps
+        self.saxpy_conserved_quantities(0.25, 0.75)
+ 
         #update ghosts
         self.update_ghosts()
@@ -1243 +1310 @@
         #Update boundary values
         self.update_boundary()
+
+        #set substep time
+        self.time = initial_time + self.timestep*0.5
+
+        #------------------------------------
+        #Third Euler step
+        #------------------------------------
+            
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine final and initial values
+        #of conserved quantities and cleanup
+        #------------------------------------
+        #combine steps
+        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
+ 
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+        #set substep time
+        self.time = initial_time + self.timestep       
+        
+
 
 
@@ -1344 +1445 @@
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
-        timestep = min(self.CFL*self.timestep, max_timestep)
+        timestep = min(self.CFL*self.flux_timestep, max_timestep)
 
         #Record maximal and minimal values of timestep for reporting

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -824 +824 @@
 
 
-    domain.timestep = timestep
+    domain.flux_timestep = timestep
 
 #MH090605 The following method belongs to the shallow_water domain class
@@ -875 +875 @@
          compute_fluxes_ext_central as compute_fluxes_ext
 
-    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon,
+    domain.flux_timestep = compute_fluxes_ext(timestep, domain.epsilon,
                                          domain.H0,
                                          domain.g,

anuga_work/development/demos/netherlands.py

@@ -83 +83 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 50
+N = 100
 
 
@@ -97 +97 @@
 domain.check_integrity()
 
+domain.set_timestepping_method('rk3')
+domain.set_default_order(2)
+
 #Setup order and all the beta's for the limiters (these should become defaults
-domain.default_order = 2
-domain.beta_w      = 1.0
-domain.beta_w_dry  = 0.2
-domain.beta_uh     = 1.0
-domain.beta_uh_dry = 0.2
-domain.beta_vh     = 1.0
-domain.beta_vh_dry = 0.2
 
-domain.alpha_balance = 100.0
+## domain.beta_w      = 1.0
+## domain.beta_w_dry  = 0.2
+## domain.beta_uh     = 1.0
+## domain.beta_uh_dry = 0.2
+## domain.beta_vh     = 1.0
+## domain.beta_vh_dry = 0.2
+
+## domain.alpha_balance = 100.0
 
 #Output params
@@ -117 +120 @@
 #Set bed-slope and friction
 inflow_stage = 0.5
-manning = 0.03
+manning = 0.0
 Z = Weir(inflow_stage)
 
@@ -186 +189 @@
 
 
-v.join()
+vis.join()