Changeset 8547 for trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

Timestamp:

Aug 31, 2012, 8:38:25 PM (13 years ago)

Author:

davies

Message:

Major experimental changes to balanced_dev

File:

• trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py (modified) (6 diffs)

trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

@@ -15 +15 @@
 #from balanced_basic import *
 from balanced_dev import *
+#from anuga_tsunami import *
 #---------
 #Setup computational domain
@@ -25 +26 @@
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
 #domain.set_store_vertices_uniquely(True)
-domain.minimum_allowed_height=0.001
-
 #------------------
 # Define topography
@@ -38 +37 @@
 
 def stagefun(x,y):
-    stge=-0.2*scale_me #-0.1*(x>0.5) -0.2*(x<=0.5)
+    stge=-0.2*scale_me # +0.01*(x>0.9) 
     #topo=topography(x,y) 
     return stge#*(stge>topo) + (topo)*(stge<=topo)
@@ -44 +43 @@
 domain.set_quantity('elevation',topography)     # Use function for elevation
 domain.get_quantity('elevation').smooth_vertex_values() 
-
-domain.set_quantity('friction',0.00)             # Constant friction
+domain.set_quantity('friction',0.03)             # Constant friction
 
 #def frict_change(x,y):
@@ -74 +72 @@
 # Associate boundary tags with boundary objects
 #----------------------------------------------
-domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
 
 #------------------------------
@@ -80 +78 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.1,finaltime=20.0):
+for t in domain.evolve(yieldstep=0.2,finaltime=20.0):
     print domain.timestepping_statistics()
     xx = domain.quantities['xmomentum'].centroid_values
     yy = domain.quantities['ymomentum'].centroid_values
     dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 
-    dd = (dd)*(dd>1.0e-06)+1.0e-03
-    vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
+    dd_raw=1.0*dd
+    dd = (dd)*(dd>1.0e-03)+1.0e-03
+    vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
     vv = vv*(dd>1.0e-03)
     print 'Peak velocity is: ', vv.max(), vv.argmax()
+    print 'Volume is', sum(dd_raw*domain.areas)    
+
+
+vel_final_inds=(vv>1.0e-01).nonzero()
+print 'vel_final_inds', vel_final_inds
 
 print 'Finished'