Changeset 5153 for anuga_work/development/demos/netherlands.py

Timestamp:

Mar 11, 2008, 4:21:13 AM (17 years ago)

Author:

ole

Message:

Updated the old netherlands example in preparation for Steve's new deep-water limiter.

File:

• anuga_work/development/demos/netherlands.py (modified) (7 diffs)

anuga_work/development/demos/netherlands.py

@@ -6 +6 @@
 """
 
-######################
-# Module imports
-#
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary, Dirichlet_boundary
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+import os
+
+#from anuga.visualiser import RealtimeVisualiser
 #import rpdb
 #rpdb.set_active()
 
-from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
-     Transmissive_boundary, Constant_height, Constant_stage
 
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from Numeric import array
-from anuga.visualiser import RealtimeVisualiser
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+
+N = 150 # size = 45000
+N = 130 # size = 33800
+N = 600 # Size = 720000
+N = 100
+
+points, elements, boundary = rectangular_cross(N, N)
+domain = Domain(points, elements, boundary, use_inscribed_circle=True)
+
+domain.check_integrity()
+
+domain.set_name(os.path.splitext(__file__)[0])
+domain.set_timestepping_method('rk3')
+domain.set_default_order(2)
+domain.set_store_vertices_uniquely(True) # Store as internally represented
+domain.tight_slope_limiters = True
+print domain.statistics()
+
+
+# Setup order and all the beta's for the limiters (these should become defaults
+
+#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
+
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
 
 class Weir:
@@ -35 +75 @@
         z = zeros(N, Float)
         for i in range(N):
-            z[i] = -x[i]/20  #General slope
+            z[i] = -x[i]/20  # General slope
 
-            #Flattish bit to the left
+            # Flattish bit to the left
             if x[i] <= 0.3:
                 #z[i] = -x[i]/5
@@ -43 +83 @@
 
 
-            #Weir
+            # Weir
             if x[i] > 0.3 and x[i] < 0.4:
                 z[i] = -x[i]/20+1.2
 
-            #Dip
+            # Dip
             #if x[i] > 0.6 and x[i] < 0.9:
             #    z[i] = -x[i]/20-0.5  #-y[i]/5
 
-            #Hole in weir
+            # Hole in weir
             #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
             if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6:
@@ -57 +97 @@
                 z[i] = -x[i]/20
 
-            #Poles
+            # Poles
             #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
             #       x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
@@ -67 +107 @@
 
 
-            #Wall
+            # Wall
             if x[i] > 0.995:
                 z[i] = -x[i]/20+0.3
@@ -74 +114 @@
 
 
+inflow_stage = 0.5
+manning = 0.0
 
-######################
-# Domain
-#
+domain.set_quantity('elevation', Weir(inflow_stage))
+domain.set_quantity('friction', manning)
+domain.set_quantity('stage', expression='elevation + 0.0')
 
 
-N = 150 #size = 45000
-N = 130 #size = 33800
-N = 600 #Size = 720000
-N = 100
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0]) # Constant inflow
+domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
 
 
-#N = 15
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
 
-print 'Creating domain'
-#Create basic mesh
-points, elements, boundary = rectangular_cross(N, N)
 
-#Create shallow water domain
-domain = Domain(points, elements, boundary, use_inscribed_circle=True)
+if N <= 150:
+    # Initialise real-time visualiser 
 
-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.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
-domain.smooth = False
-domain.reduction = min  #Looks a lot better on top of steep slopes
-
-print "Number of triangles = ", len(domain)
-print "Extent = ", domain.get_extent()
-
-#Set bed-slope and friction
-inflow_stage = 0.5
-manning = 0.0
-Z = Weir(inflow_stage)
-
-if N > 150:
-    domain.checkpoint = False
-    domain.store = True    #Store for visualisation purposes
-    domain.format = 'sww'   #Native netcdf visualisation format
-    import sys, os
-    #FIXME: This was os.path.splitext but caused weird filenames based on root
-    base = os.path.basename(sys.argv[0])
-    basename, _ = os.path.splitext(base)
-    domain.set_name(basename)
-else:
-    domain.checkpoint = False
-    domain.store = False
-    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
+    #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()
 
 
 
-print 'Field values'
-domain.set_quantity('elevation', Z)
-domain.set_quantity('friction', manning)
-
-
-######################
-# Boundary conditions
-#
-print 'Boundaries'
-Br = Reflective_boundary(domain)
-Bt = Transmissive_boundary(domain)
-
-#Constant inflow
-Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
-
-
-#Set boundary conditions
-domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
-
-
-######################
-#Initial condition
-#
-print 'Initial condition'
-domain.set_quantity('stage', expression='elevation + 0.0')
-
-#Evolve
 import time
 t0 = time.time()
 
 for t in domain.evolve(yieldstep = 0.005, finaltime = None):
-    domain.write_time()
-    #domain.write_boundary()
-    vis.update()
+    print domain.timestepping_statistics()
     print domain.quantities['stage'].get_values(location='centroids',
                                                 indices=[0])
+                                                
+    #vis.update()                                               
     #time.sleep(0.1)
     #raw_input('pause>')
@@ -184 +162 @@
     #rpdb.set_active()
     #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
-vis.evolveFinished()
+#vis.evolveFinished()
 
 print 'That took %.2f seconds' %(time.time()-t0)
 
-
-vis.join()
+#vis.join()