Changeset 5741

Timestamp:

Sep 5, 2008, 9:16:45 PM (17 years ago)

Author:

steve

Message:

Updated timestepping

Location:

anuga_work/development

Files:

• anuga_1d/config.py (modified) (1 diff)
• anuga_1d/dam_h_elevation.py (modified) (2 diffs)
• anuga_1d/domain.py (modified) (6 diffs)
• anuga_1d/quantity.py (modified) (2 diffs)
• anuga_1d/shallow_water_domain.py (modified) (3 diffs)
• anuga_1d/shallow_water_domain_new.py (modified) (6 diffs)
• anuga_1d/test_shallow_water.py (modified) (5 diffs)
• demos/netherlands.py (modified) (5 diffs)
• shallow_water_1d_old/quantity.py (modified) (3 diffs)

anuga_work/development/anuga_1d/config.py

@@ -11 +11 @@
 
 min_timestep = 1.0e-6 #Should be computed based on geometry
-max_timestep = 1000
+max_timestep = 1.0e3
+max_smallsteps = 50  # Max number of degenerate steps allowed b4 trying first order
 #This is how:
 #Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)

anuga_work/development/anuga_1d/dam_h_elevation.py

@@ -1 +1 @@
 import os
 from math import sqrt, sin, cos, pi, exp
-from shallow_water_domain import *
+from shallow_water_domain_new import *
 from Numeric import zeros, Float
 #from analytic_dam_sudi import AnalyticDam
@@ -42 +42 @@
 
 domain.default_order = 2
-domain.default_time_order = 1
+domain.set_timestepping_method('euler')
+#domain.default_time_order = 2
 domain.cfl = 1.0
 domain.limiter = "vanleer"

anuga_work/development/anuga_1d/domain.py

@@ -724 +724 @@
     #--------------------------    
 
-    def evolve(self, yieldstep = None, finaltime = None, duration = None,
+    def evolve(self, yieldstep = None,
+               finaltime = None,
+               duration = None,
                skip_initial_step = False):
         """Evolve model through time starting from self.starttime.
@@ -759 +761 @@
 
         from config import min_timestep, max_timestep, epsilon
-
-        print 'yieldstep', yieldstep
-        print 'duration', duration
 
         # FIXME: Maybe lump into a larger check prior to evolving
@@ -1388 +1387 @@
     def update_timestep(self, yieldstep, finaltime):
 
-        from config import min_timestep
+        from config import min_timestep, max_timestep
 
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
-        timestep = self.CFL*self.timestep
+        timestep = min(self.CFL*self.flux_timestep, max_timestep)
 
         #Record maximal and minimal values of timestep for reporting
@@ -1459 +1458 @@
     
     #def update_conserved_quantities(self):
-    def update_conserved_quantities(self,timestep):
+    def update_conserved_quantities(self):
         """Update vectors of conserved quantities using previously
         computed fluxes specified forcing functions.
@@ -1469 +1468 @@
         d = len(self.conserved_quantities)
 
-        #timestep = self.timestep
+        timestep = self.timestep
 
         #Compute forcing terms
-        #self.compute_forcing_terms()
+        self.compute_forcing_terms()
 
         #Update conserved_quantities
@@ -1479 +1478 @@
             Q.update(timestep)
 
-            #Clean up
-            #Note that Q.explicit_update is reset by compute_fluxes
-
-            #MH090605 commented out the following since semi_implicit_update is now re-initialized
-            #at the end of the _update function in quantity_ext.c (This is called by the
-            #preceeding Q.update(timestep) statement above).
-            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
-            #to 8.35 secs
-
-            #Q.semi_implicit_update[:] = 0.0
+
 
 if __name__ == "__main__":

anuga_work/development/anuga_1d/quantity.py

@@ -861 +861 @@
         #backup_centroid_values(self)
         N = self.domain.number_of_elements
-        for i in range(self.N):
-            quantity.centroid_backup_values[i] = quantity.centroid_values[i]
+        for i in range(N):
+            self.centroid_backup_values[i] = self.centroid_values[i]
 
     def saxpy_centroid_values(self,a,b):
@@ -869 +869 @@
         #saxpy_centroid_values(self,a,b)
         N = self.domain.number_of_elements
-        for i in range(self.N):
-            quantity.centroid_backup_values[i] = a*quantity.centroid_values[i] + b*quantity.centroid_backup_values[i]        
+        for i in range(N):
+            self.centroid_backup_values[i] = a*self.centroid_values[i] + b*self.centroid_backup_values[i]        
 
 class Conserved_quantity(Quantity):

anuga_work/development/anuga_1d/shallow_water_domain.py

@@ -575 +575 @@
         #print "flux cell",k,flux[0]
 
-    domain.timestep = timestep
+    domain.flux_timestep = timestep
     #print domain.quantities['stage'].centroid_values
 
@@ -652 +652 @@
         from comp_flux_ext import compute_fluxes_ext
         
-        domain.timestep = compute_fluxes_ext(timestep,
+        domain.flux_timestep = compute_fluxes_ext(timestep,
                                   epsilon,
                                   g,
@@ -685 +685 @@
         from comp_flux_ext import compute_fluxes_ext_short
 
-        domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
 
 

anuga_work/development/anuga_1d/shallow_water_domain_new.py

@@ -165 +165 @@
         distribute_to_vertices_and_edges(self)
         
-    def evolve(self, yieldstep = None, finaltime = None,
+    def evolve(self, yieldstep = None,
+               finaltime = None,
+               duration = None,
                skip_initial_step = False):
         """Specialisation of basic evolve method from parent class
@@ -172 +174 @@
         self.distribute_to_vertices_and_edges()
         
-        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
+        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration,
                                        skip_initial_step):
             #Pass control on to outer loop for more specific actions
@@ -439 +441 @@
         
 
-    domain.timestep = timestep
+    domain.flux_timestep = timestep
     #print domain.quantities['stage'].centroid_values
 
@@ -468 +470 @@
         from comp_flux_ext import compute_fluxes_ext
         
-        domain.timestep = compute_fluxes_ext(timestep,
+        domain.flux_timestep = compute_fluxes_ext(timestep,
                                   epsilon,
                                   g,
@@ -501 +503 @@
         from comp_flux_ext import compute_fluxes_ext_short
 
-        domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
 
 
@@ -526 +528 @@
         max_speed_array = domain.max_speed_array
                 
-        from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced                  #from comp_flux_ext import compute_fluxes_ext
+        from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
         
-        domain.timestep = compute_fluxes_ext_wellbalanced(timestep,
+        domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
                                   epsilon,
                                   g,

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -5 +5 @@
 
 
-from shallow_water_domain import *
-from shallow_water_domain import flux_function as domain_flux_function
+from shallow_water_domain_new import *
+from shallow_water_domain_new import flux_function as domain_flux_function
 
 from Numeric import allclose, array, ones, Float, maximum, zeros
@@ -41 +41 @@
         #print domain.quantities['stage'].explicit_update
         #print domain.quantities['xmomentum'].explicit_update
-        #print stage_ud
+        print stage_ud
 
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
@@ -133 +133 @@
 
         domain.default_order = 1
-        domain.default_time_order = 1
+        domain.set_timestepping_method('euler')
+        
         yieldstep=1.0
         finaltime=1.0
@@ -167 +168 @@
 
         domain.default_order = 2
-        domain.default_time_order = 1
+        domain.set_timestepping_method('rk2')
         yieldstep=1.0
         finaltime=1.0
@@ -191 +192 @@
 
         domain.default_order = 2
-        domain.default_time_order = 2
+        domain.set_timestepping_method('rk3')
+
         yieldstep=1.0
         finaltime=1.0

anuga_work/development/demos/netherlands.py

@@ -27 +27 @@
 N = 130 # size = 33800
 N = 600 # Size = 720000
-N = 100
+N = 50
 
 points, elements, boundary = rectangular_cross(N, N)
@@ -34 +34 @@
 domain.check_integrity()
 
+import sys
+base = os.path.basename(sys.argv[0])
+domain.simulation_name = 'netherlands'
 domain.set_name(os.path.splitext(__file__)[0])
-domain.set_timestepping_method('rk3')
+domain.set_timestepping_method('rk2')
 domain.set_default_order(2)
 domain.set_store_vertices_uniquely(True) # Store as internally represented
@@ -130 +133 @@
 domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
 
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = 'dam_break_'+time
+output_file = 'dam_break'
+
+copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
 
 #------------------------------------------------------------------------------
@@ -139 +154 @@
     # Initialise real-time visualiser 
 
-    pass
-    #vis = RealtimeVisualiser(domain)
-    #vis.render_quantity_height("elevation", dynamic=False)
-    #vis.render_quantity_height("stage", dynamic=True)
-    #vis.colour_height_quantity('stage', (0.0, 0.0, 0.8))
-    #vis.start()
+    #pass
+    visualise = True
+    if visualise:
+        from anuga.visualiser import RealtimeVisualiser
+        vis = RealtimeVisualiser(domain)
+        vis.render_quantity_height("elevation", offset=0.01, dynamic=False)
+        vis.render_quantity_height("stage", dynamic=True)
+        vis.colour_height_quantity('stage', (0.2, 0.2, 0.8))
+        vis.start()
+        import time
+        time.sleep(2.0)
+    
+    
 
 
@@ -151 +173 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.005, finaltime = None):
+for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
     print domain.timestepping_statistics()
-    print domain.quantities['stage'].get_values(location='centroids',
-                                                indices=[0])
-                                                
-    #vis.update()                                               
-    #time.sleep(0.1)
-    #raw_input('pause>')
-    #V.update_quantity('stage')
-    #rpdb.set_active()
-    #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
-#vis.evolveFinished()
+    #print domain.quantities['stage'].get_values(location='centroids',
+                                                #indices=[0])
+    if visualise:
+        vis.update()                                            
+
+        
+
+if visualise: vis.evolveFinished()
 
 print 'That took %.2f seconds' %(time.time()-t0)

anuga_work/development/shallow_water_1d_old/quantity.py

@@ -56 +56 @@
         #Intialise centroid values
         self.interpolate()
+
+
+        from Numeric import zeros, Float
+        
+        #Allocate space for boundary values
+        #L = len(domain.boundary)
+        self.boundary_values = zeros(2, Float) #assumes no parrellism
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+        
+        N = domain.number_of_elements
+        self.explicit_update = zeros(N, Float )
+        self.semi_implicit_update = zeros(N, Float )
+
+        self.gradients = zeros(N, Float)
+        self.qmax = zeros(self.centroid_values.shape, Float)
+        self.qmin = zeros(self.centroid_values.shape, Float)
+
+        self.beta = domain.beta        
 
     #Methods for operator overloading
@@ -381 +401 @@
         return integral
 
-class Conserved_quantity(Quantity):
-    """Class conserved quantity adds to Quantity:
-
-    storage and method for updating, and
-    methods for extrapolation from centropid to vertices inluding
-    gradients and limiters
-    """
-
-    def __init__(self, domain, vertex_values=None):
-        Quantity.__init__(self, domain, vertex_values)
-
-        from Numeric import zeros, Float
-        
-        #Allocate space for boundary values
-        #L = len(domain.boundary)
-        self.boundary_values = zeros(2, Float) #assumes no parrellism
-
-        #Allocate space for updates of conserved quantities by
-        #flux calculations and forcing functions
-        
-        N = domain.number_of_elements
-        self.explicit_update = zeros(N, Float )
-        self.semi_implicit_update = zeros(N, Float )
-
-        self.gradients = zeros(N, Float)
-        self.qmax = zeros(self.centroid_values.shape, Float)
-        self.qmin = zeros(self.centroid_values.shape, Float)
-
-        self.beta = domain.beta
 
     def update(self, timestep):
@@ -863 +854 @@
                 else:
                     qv[k,i] = qc[k]
+    
+
+class Conserved_quantity(Quantity):
+    """Class conserved quantity adds to Quantity:
+
+    storage and method for updating, and
+    methods for extrapolation from centropid to vertices inluding
+    gradients and limiters
+    """
+
+    def __init__(self, domain, vertex_values=None):
+        Quantity.__init__(self, domain, vertex_values)
+
+        print "Use Quantity instead of Conserved_quantity"
+