Changeset 3697

Timestamp:

Oct 5, 2006, 12:21:52 PM (17 years ago)

Author:

ole

Message:

Played with limiter example: run_sw_rectangle.py

File:

• anuga_core/source/anuga/examples/run_sw_rectangle.py (modified) (4 diffs)

anuga_core/source/anuga/examples/run_sw_rectangle.py

@@ -1 +1 @@
 #!/usr/bin/env python
-#########################################################
-#
-#  Example showing improved wet dry limiting.
-#
-#
-#  Authors: Steve Roberts
-# October
-#
-#
-#
-#########################################################
+"""
+  Example showing improved wet dry limiting.
+  This script runs the same simulation twice and stores the outputs in
+  old_limiter.sww and new_limiter.sww respectively.
+  
+  Authors: Steve Roberts
+  October
+"""
 
-
+#-----------------------------------------------------------------
+# Common structures
+#-----------------------------------------------------------------
+import time
 from Numeric import array
 from anuga.shallow_water import Domain
-
-
-# mesh partition routines
+from anuga.shallow_water import Reflective_boundary
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
-
-
-M = 30
-points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0)
-domain = Domain(points, vertices, boundary)
-print 'number of triangles = ', domain.number_of_elements
-
-
-#---------------------------
-#Boundaries
-#---------------------------
-from anuga.shallow_water import Reflective_boundary
-
-R = Reflective_boundary(domain)
-
-
-domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
-domain.check_integrity()
 
 class Set_Stage:
@@ -53 +33 @@
         return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
 
+M = 30
+points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0)
+
+yieldstep = 0.002
+finaltime = 0.06
+rect = [0.0, 0.0, 1.0, 1.0]
+
+
+#-----------------------------------------------------------------
+# Create domain for "old limiter" scenario
+#-----------------------------------------------------------------
+domain = Domain(points, vertices, boundary)
+domain.set_name('old_limiter_second_order')
+print 'Number of triangles =', len(domain)
+
+# Turn on the visualisation
+try:
+    domain.initialise_visualiser()
+except:
+    pass
+
+
+#-----------------------------------------------------------------
+# Boundaries and Initial conditions
+#-----------------------------------------------------------------
+R = Reflective_boundary(domain)
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
 
-import time
-t0 = time.time()
-
-
-# Turn on the visualisation
-
-rect = [0.0, 0.0, 1.0, 1.0]
-domain.initialise_visualiser()
-
-yieldstep = 0.002
-finaltime = 0.05
-
-
-#===============================================================================
-#Old Limiter
-#===============================================================================
-
+#-----------------------------------------------------------------
+# Values for old limiter
+#-----------------------------------------------------------------
 domain.default_order = 2
 domain.beta_w      = 0.9
@@ -80 +73 @@
 domain.beta_vh_dry = 0.9
 
-#Check that the boundary value gets propagated to all elements
+#-----------------------------------------------------------------
+# Evolve
+#-----------------------------------------------------------------
+t0 = time.time()
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
 
-
 print 'That took %.2f seconds' %(time.time()-t0)
 print 'Note the small timesteps and the irregular flow'
-raw_input('press return to continue')
+#raw_input('press return to continue')
 
-#===============================================================================
-#New Limiter
-#===============================================================================
 
-t0 = time.time()
+#-----------------------------------------------------------------
+# Create domain for "new limiter" scenario (2 order)
+#-----------------------------------------------------------------
+domain = Domain(points, vertices, boundary)
+domain.set_name('new_limiter_second_order')
+print 'Number of triangles =', len(domain)
+
+# Turn on the visualisation
+try:
+    domain.initialise_visualiser()
+except:
+    pass
+
+#-----------------------------------------------------------------
+# Boundaries and Initial conditions
+#-----------------------------------------------------------------
+R = Reflective_boundary(domain)
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
-domain.set_quantity('xmomentum', 0.0)
-domain.set_quantity('ymomentum', 0.0)
-domain.time = 0.0
-domain.default_order = 2
+
+#-----------------------------------------------------------------
+# Values for new limiter
+#-----------------------------------------------------------------
+domain.set_default_order(2)
 domain.minimum_allowed_height = 0.001
 domain.beta_w      = 1.0
@@ -107 +117 @@
 domain.beta_vh_dry = 0.2
 
-#Check that the boundary value gets propagated to all elements
+#-----------------------------------------------------------------
+# Evolve
+#-----------------------------------------------------------------
+t0 = time.time()
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
 
-
 print 'That took %.2f seconds' %(time.time()-t0)
 print 'Note more uniform and large timesteps'
-raw_input('press return to continue')
-
-#===============================================================================
-#First Order
-#===============================================================================
-
-t0 = time.time()
-domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
-domain.set_quantity('xmomentum', 0.0)
-domain.set_quantity('ymomentum', 0.0)
-domain.time = 0.0
-
-domain.default_order = 1
+#raw_input('press return to continue')
 
 
-#Check that the boundary value gets propagated to all elements
+#-----------------------------------------------------------------
+# Create domain for "new limiter" scenario (1 order)
+#-----------------------------------------------------------------
+domain = Domain(points, vertices, boundary)
+domain.set_name('new_limiter_first_order')
+print 'Number of triangles =', len(domain)
+
+# Turn on the visualisation
+try:
+    domain.initialise_visualiser()
+except:
+    pass
+
+#-----------------------------------------------------------------
+# Boundaries and Initial conditions
+#-----------------------------------------------------------------
+R = Reflective_boundary(domain)
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
+domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+
+#-----------------------------------------------------------------
+# Values for new limiter first order
+#-----------------------------------------------------------------
+domain.set_default_order(1)
+domain.minimum_allowed_height = 0.001
+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
+
+#-----------------------------------------------------------------
+# Evolve
+#-----------------------------------------------------------------
+t0 = time.time()
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
 
-
 print 'That took %.2f seconds' %(time.time()-t0)
 
+