Changeset 8204 for trunk/anuga_work/development

Timestamp:

Sep 1, 2011, 7:28:58 PM (14 years ago)

Author:

steve

Message:

Getting the 1d tsunami model working

Location:

trunk/anuga_work/development/2010-projects/anuga_1d

Files:

• channel/channel_flow_2_padarn.py (modified) (6 diffs)
• sww/run_dry_dam.py (modified) (1 diff)
• sww/sww_domain.py (modified) (2 diffs)

trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_2_padarn.py

@@ -1 +1 @@
 import os
 from math import sqrt, pow, pi 
-from channel_domain import *
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
-from config import g, epsilon
+
+import numpy as np
+
+from anuga_1d.channel.channel_domain import *
+from anuga_1d.config import g, epsilon
+from anuga_1d.base.generic_mesh import uniform_mesh
 
 
@@ -19 +22 @@
 
 def bed(x):
-    y = zeros(len(x),Float)
+    y = np.zeros(len(x),np.float)
     for i in range(len(x)):
         if x[i]<525 and x[i]>475:
@@ -46 +49 @@
 print "Evaluating domain with %d cells" %N
 cell_len = L/N # Origin = 0.0
-points = zeros(N+1,Float)
+points = np.zeros(N+1,np.float)
 
 # Define the centroid points
@@ -62 +65 @@
 
 # Set boundry type, order, timestepping method and limiter
-domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+#domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+domain.set_boundary({'exterior': Reflective_boundary(domain)})
 domain.order = 2
 domain.set_timestepping_method('rk2')
@@ -72 +76 @@
 t0 = time.time()
 
+finaltime= 0.0
+
+#===================================================================
+# Time loop
+#===================================================================
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
@@ -84 +93 @@
 VelocityQ = domain.quantities['velocity'].vertex_values
 x = X.flat
-z = domain.quantities['elevation'].vertex_values.flat
-stage=HeightQ.flat+z
+z = domain.quantities['elevation'].vertex_values
+stage=HeightQ+z
+stage = stage.flat
+z = z.flat
 
 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot

trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

@@ -81 +81 @@
 k = 0
 
-N = 10000
+N = 2000
 print "Evaluating domain with %d cells" %N
 domain = Domain(*uniform_mesh(N,x_1=L))

trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -235 +235 @@
         
     #Calculate height (and fix negatives)better be non-negative!
+    w_C[:] = numpy.maximum(w_C, z_C)
     h_C[:] = w_C - z_C
+
+
     u_C[:]  = uh_C/(h_C + h0/h_C)
 
@@ -282 +285 @@
     #print 'any' ,numpy.any( h_V[:,1] < 0.0)
 
-    h00 = 1e-12
+    #h00 = 1e-12
     #print h00
 
-    h_V[:,:] = numpy.where (h_V <= h00, 0.0, h_V)
-    u_V[:,:] = numpy.where (h_V <= h00, 0.0, u_V)
+    #h_V[:,:] = numpy.where (h_V <= h00, 0.0, h_V)
+    #u_V[:,:] = numpy.where (h_V <= h00, 0.0, u_V)
 
 #    # protect from edge values going negative