Changeset 7842

Timestamp:

Jun 15, 2010, 5:34:24 PM (14 years ago)

Author:

steve

Message:

Adding in a few new files

Location:

anuga_work/development/2010-projects/anuga_1d

Files:

• channel/channel_domain.py (modified) (26 diffs)
• channel/channel_domain_ext.c (modified) (1 diff)
• channel/channel_flow_1_padarn.py (modified) (7 diffs)
• channel/profile_channel.py (added)
• channel/test_padarn.py (modified) (8 diffs)
• generic/generic_domain.py (modified) (3 diffs)
• generic/generic_mesh.py (added)

anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

@@ -3 +3 @@
 the shallow water wave equation.
 
-This module contains a specialisation of class Domain from module domain.py
-consisting of methods specific to the Shallow Water Wave Equation
+This module contains a specialisation of class Generic_domain from module
+generic_domain.py
+consisting of methods specific to Channel flow using the Shallow Water Wave Equation
 
 This particular modification of the Domain class implements the ability to
@@ -14 +15 @@
 
 where
-------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
-U = [A, Q]
-E = [Q, Q^2/A + gh^2/2]
+
+U = [A, Q] = [b*h, u*b*h]
+E = [Q, Q^2/A + g*b*h^2/2]
 S represents source terms forcing the system
-(e.g. gravity, friction, wind stress, ...)
+(e.g. gravity, boundary_stree, friction, wind stress, ...)
+gravity = -g*b*h*z_x
+boundary_stress = 1/2*g*b_x*h^2
 
 and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
@@ -25 +28 @@
 
 symbol    variable name    explanation
+A         area             Wetted area = b*h
+Q         discharge        flux of water = u*b*h
 x         x                horizontal distance from origin [m]
 z         elevation        elevation of bed on which flow is modelled [m]
@@ -31 +36 @@
 u                          speed in the x direction [m/s]
 uh        xmomentum        momentum in the x direction [m^2/s]
-
+b         width            width of channel
 eta                        mannings friction coefficient [to appear]
 nu                         wind stress coefficient [to appear]
 
-The conserved quantities are w, uh
+The conserved quantities are A, Q
 --------------------------------------------------------------------------
 For details see e.g.
@@ -48 +53 @@
 
 
-from domain import *
-Generic_Domain = Domain #Rename
+from anuga_1d.generic.generic_domain import *
+import numpy
+
 
 #Shallow water domain
-class Domain(Generic_Domain):
+class Domain(Generic_domain):
 
     def __init__(self, coordinates, boundary = None, tagged_elements = None):
@@ -59 +65 @@
         evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
         other_quantities = ['friction']
-        Generic_Domain.__init__(self,
+        Generic_domain.__init__(self,
                                 coordinates = coordinates,
                                 boundary = boundary,
@@ -67 +73 @@
                                 tagged_elements = tagged_elements)
         
-        from config import minimum_allowed_height, g, h0
+        from anuga_1d.config import minimum_allowed_height, g, h0
         self.minimum_allowed_height = minimum_allowed_height
         self.g = g
@@ -74 +80 @@
         self.discontinousb = False
 
-     
+
+        # forcing terms gravity and boundary stress are included in the flux calculation
         #self.forcing_terms.append(gravity)
         #self.forcing_terms.append(boundary_stress)
@@ -88 +95 @@
         
         #Reduction operation for get_vertex_values
-        from util import mean
+        from anuga_1d.generic.util import mean
         self.reduction = mean
         #self.reduction = min  #Looks better near steep slopes
@@ -94 +101 @@
         self.set_quantities_to_be_stored(['area','discharge'])
 
-        self.__doc__ = 'channel_domain_Ab'
+        self.__doc__ = 'channel_domain'
 
         self.check_integrity()
@@ -118 +125 @@
         msg = 'Fifth evolved quantity must be "velocity"'
         assert self.evolved_quantities[4] == 'velocity', msg
-        msg = 'Fifth evolved quantity must be "width"'
+        msg = 'Sixth evolved quantity must be "width"'
         assert self.evolved_quantities[5] == 'width', msg
-
-        Generic_Domain.check_integrity(self)
+        msg = 'Seventh evolved quantity must be "stage"'
+        assert self.evolved_quantities[6] == 'stage', msg
+
+        Generic_domain.check_integrity(self)
 
     def compute_fluxes(self):
@@ -140 +149 @@
 #-----------------------------------
 def compute_fluxes_channel(domain):
-    from Numeric import zeros, Float
     import sys
-
-        
     timestep = float(sys.maxint)
 
-    area    = domain.quantities['area']
-    discharge     = domain.quantities['discharge']
-    bed      = domain.quantities['elevation']
-    height   = domain.quantities['height']
-    velocity = domain.quantities['velocity']
-    width    = domain.quantities['width']
-
-
-    from channel_domain_ext import compute_fluxes_channel_ext
+    area       = domain.quantities['area']
+    discharge  = domain.quantities['discharge']
+    bed        = domain.quantities['elevation']
+    height     = domain.quantities['height']
+    velocity   = domain.quantities['velocity']
+    width      = domain.quantities['width']
+
+
+    from anuga_1d.channel.channel_domain_ext import compute_fluxes_channel_ext
     domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
 
@@ -168 +174 @@
 
     import sys
-    from Numeric import zeros, Float
-    from config import epsilon, h0
-   ##  linearb(domain)
-
-
+    from anuga_1d.config import epsilon, h0
 
     
@@ -178 +180 @@
 
     #Shortcuts
-    Area = domain.quantities['area']
-    Discharge = domain.quantities['discharge']
-    Bed = domain.quantities['elevation']
-    Height = domain.quantities['height']
-    Velocity = domain.quantities['velocity']
-    Width = domain.quantities['width']
-    Stage = domain.quantities['stage']
+    area      = domain.quantities['area']
+    discharge = domain.quantities['discharge']
+    bed       = domain.quantities['elevation']
+    height    = domain.quantities['height']
+    velocity  = domain.quantities['velocity']
+    width     = domain.quantities['width']
+    stage     = domain.quantities['stage']
 
     #Arrays   
-    a_C   = Area.centroid_values    
-    d_C  = Discharge.centroid_values    
-    z_C   = Bed.centroid_values
-    h_C   = Height.centroid_values
-    u_C   = Velocity.centroid_values
-    b_C   = Width.centroid_values
-    w_C   = Stage.centroid_values
+    a_C   = area.centroid_values
+    d_C   = discharge.centroid_values
+    z_C   = bed.centroid_values
+    h_C   = height.centroid_values
+    u_C   = velocity.centroid_values
+    b_C   = width.centroid_values
+    w_C   = stage.centroid_values
 
     if domain.setstageflag:
-        for i in range(len(a_C)):
-            a_C[i]=(w_C[i]-z_C[i])*b_C[i]
-     
+        a_C[:,] = (w_C-z_C)*b_C
         domain.setstageflag = False
 
@@ -205 +205 @@
         
 
-    h0 = 1.0e-12        
-    #print id(h_C)
-    for i in range(N):
-                                                        
-        if a_C[i] <= h0:
-            a_C[i] = 0.0
-            h_C[i] = 0.0
-            d_C[i] = 0.0
-            u_C[i] = 0.0
-            w_C[i] = z_C[i]
-           
-            
-
-        else:
-
-            if b_C[i]<=h0:
-                a_C[i] = 0.0
-                h_C[i] = 0.0
-                d_C[i] = 0.0
-                u_C[i] = 0.0
-                w_C[i] = z_C[i]
-
-            else:
-                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
-                w_C[i] = h_C[i]+z_C[i]
-                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
+    # FIXME (SGR): Replace with numpy.where to speed up code
+    h0 = 1.0e-12
+
+    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
+    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
+
+    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
+    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
+    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
+
+    w_C[:] = h_C + z_C
+    
+#    for i in range(N):
+#
+#        if a_C[i] <= h0:
+#           a_C[i] = 0.0
+#           h_C[i] = 0.0
+#            d_C[i] = 0.0
+#            u_C[i] = 0.0
+#            w_C[i] = z_C[i]
+#        else:
+#
+#            if b_C[i]<=h0:
+#                a_C[i] = 0.0
+#                h_C[i] = 0.0
+#                d_C[i] = 0.0
+#                u_C[i] = 0.0
+#                w_C[i] = z_C[i]
+#
+#            else:
+#                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
+#                w_C[i] = h_C[i]+z_C[i]
+#                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
                 
-   
-   
-        
+
+    bed.extrapolate_second_order()
+    
     for name in ['velocity','stage']:
         Q = domain.quantities[name]
@@ -243 +250 @@
         else:
             raise 'Unknown order'
-    a_V  = domain.quantities['area'].vertex_values
-    w_V  = domain.quantities['stage'].vertex_values                     
-    z_V  = domain.quantities['elevation'].vertex_values 
-    h_V  = domain.quantities['height'].vertex_values
-    u_V  = domain.quantities['velocity'].vertex_values          
-    d_V = domain.quantities['discharge'].vertex_values
-    b_V = domain.quantities['width'].vertex_values
+
+
+    a_V  = area.vertex_values
+    w_V  = stage.vertex_values
+    z_V  = bed.vertex_values
+    h_V  = height.vertex_values
+    u_V  = velocity.vertex_values
+    d_V  = discharge.vertex_values
+    b_V  = width.vertex_values
 
    
-    
-    for i in range(len(h_C)):
-        for j in range(2):
-##             if b_V[i,j] < h0 :
-##                 a_V[i,j]=0
-##                 h_V[i,j]=0
-##                 d_V[i,j]=0
-##                 u_V[i,j]=0
-##             else:
-                
-            h_V[i,j] = w_V[i,j]-z_V[i,j]
-            if h_V[i,j]<h0:
-                h_V[i,j]=0
-                w_V[i,j]=z_V[i,j]
-            a_V[i,j] = b_V[i,j]*h_V[i,j]
-            d_V[i,j]=u_V[i,j]*a_V[i,j]
-            
-                
-                
-    
-
+    #FIXME (SGR): replace with numpy.where
+
+
+    h_V[:] = w_V-z_V
+
+    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
+    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
+    a_V[:] = b_V*h_V
+    d_V[:] = u_V*a_V
+
+#    for i in range(len(h_C)):
+#        for j in range(2):
+#            h_V[i,j] = w_V[i,j]-z_V[i,j]
+#            if h_V[i,j]<h0:
+#                h_V[i,j]=0
+#                w_V[i,j]=z_V[i,j]
+#            a_V[i,j] = b_V[i,j]*h_V[i,j]
+#            d_V[i,j]=u_V[i,j]*a_V[i,j]
     
     return
@@ -278 +284 @@
 
 #--------------------------------------------------------
-#Boundaries - specific to the shallow_water_vel_domain
+#Boundaries - specific to the channel_domain
 #--------------------------------------------------------
 class Reflective_boundary(Boundary):
@@ -298 +304 @@
         #Handy shorthands
         self.normals  = domain.normals
-        self.area    = domain.quantities['area'].vertex_values
+        self.area     = domain.quantities['area'].vertex_values
         self.discharge     = domain.quantities['discharge'].vertex_values
         self.bed      = domain.quantities['elevation'].vertex_values
@@ -306 +312 @@
         self.stage    = domain.quantities['stage'].vertex_values
 
-        from Numeric import zeros, Float
-        #self.conserved_quantities = zeros(3, Float)
-        self.evolved_quantities = zeros(7, Float)
+        self.evolved_quantities = numpy.zeros(7, numpy.float)
 
     def __repr__(self):
@@ -320 +324 @@
 
         q = self.evolved_quantities
-        q[0] = self.area[vol_id, edge_id]
+        q[0] =  self.area[vol_id, edge_id]
         q[1] = -self.discharge[vol_id, edge_id]
-        q[2] = self.bed[vol_id, edge_id]
-        q[3] = self.height[vol_id, edge_id]
+        q[2] =  self.bed[vol_id, edge_id]
+        q[3] =  self.height[vol_id, edge_id]
         q[4] = -self.velocity[vol_id, edge_id]
-        q[5] = self.width[vol_id,edge_id]
-        q[6] = self.stage[vol_id,edge_id]
-
-        #print "In Reflective q ",q
-
+        q[5] =  self.width[vol_id,edge_id]
+        q[6] =  self.stage[vol_id,edge_id]
 
         return q
@@ -348 +349 @@
             raise msg
 
-        from Numeric import array, Float
-        self.evolved_quantities=array(evolved_quantities).astype(Float)
+        assert len(evolved_quantities) == 7
+
+        self.evolved_quantities=numpy.array(evolved_quantities,numpy.float)
 
     def __repr__(self):
@@ -370 +372 @@
 
 
-    Area     = domain.quantities['area']
+    Area      = domain.quantities['area']
     Discharge  = domain.quantities['discharge']
     Elevation = domain.quantities['elevation']
@@ -402 +404 @@
 def boundary_stress(domain):
 
-
     from util import gradient
     from Numeric import zeros, Float, array, sum
-
-
 
     Area     = domain.quantities['area']
@@ -415 +414 @@
 
     discharge_ud  = Discharge.explicit_update
-
-
-       
+ 
     h = Height.vertex_values
     b = Width.vertex_values
@@ -451 +448 @@
 
     uh = domain.quantities['xmomentum'].centroid_values
-    #vh = domain.quantities['ymomentum'].centroid_values
     eta = domain.quantities['friction'].centroid_values
 
     xmom_update = domain.quantities['xmomentum'].semi_implicit_update
-    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
     N = domain.number_of_elements
@@ -502 +497 @@
 
 
-def check_forcefield(f):
-    """Check that f is either
-    1: a callable object f(t,x,y), where x and y are vectors
-       and that it returns an array or a list of same length
-       as x and y
-    2: a scalar
-    """
-
-    from Numeric import ones, Float, array
-
-
-    if callable(f):
-        #N = 3
-        N = 2
-        #x = ones(3, Float)
-        #y = ones(3, Float)
-        x = ones(2, Float)
-        #y = ones(2, Float)
-        
-        try:
-            #q = f(1.0, x=x, y=y)
-            q = f(1.0, x=x)
-        except Exception, e:
-            msg = 'Function %s could not be executed:\n%s' %(f, e)
-            #FIXME: Reconsider this semantics
-            raise msg
-
-        try:
-            q = array(q).astype(Float)
-        except:
-            msg = 'Return value from vector function %s could ' %f
-            msg += 'not be converted into a Numeric array of floats.\n'
-            msg += 'Specified function should return either list or array.'
-            raise msg
-
-        #Is this really what we want?
-        msg = 'Return vector from function %s ' %f
-        msg += 'must have same lenght as input vectors'
-        assert len(q) == N, msg
-
-    else:
-        try:
-            f = float(f)
-        except:
-            msg = 'Force field %s must be either a scalar' %f
-            msg += ' or a vector function'
-            raise msg
-    return f
 
 def linearb(domain):

anuga_work/development/2010-projects/anuga_1d/channel/channel_domain_ext.c

@@ -1 +1 @@
 #include "Python.h"
-#include "Numeric/arrayobject.h"
+#include "numpy/arrayobject.h"
 #include "math.h"
 #include <stdio.h>

anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

@@ -2 +2 @@
 from math import sqrt, pow, pi
 
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+import numpy
 
-from channel_domain import *
-from config import g, epsilon
+from anuga_1d.channel.channel_domain import *
+from anuga_1d.config import g, epsilon
+from anuga_1d.generic.generic_mesh import interval_mesh
 
 
@@ -21 +22 @@
 
 def bed(x):
-    y = zeros(len(x),Float)
-    for i in range(len(x)):
-        if x[i]<525 and x[i]>475:
-            y[i]=1
-        else:
-            y[i]=0
-    return y
+    y = numpy.ones(len(x),numpy.float)
+
+    return numpy.where( (x<525) & (x>475),y,0.0)
+    
 
 def initial_discharge(x):
@@ -35 +33 @@
 
 # Set final time and yield time for simulation
-finaltime = 10.0
-yieldstep = finaltime
+finaltime = 50.0
+yieldstep = 10.0
 
 # Length of channel (m)
 L = 1000.0   
 # Define the number of cells
-number_of_cells = [200]
-
-# Define cells for finite volume and their size
-N = int(number_of_cells[0])     
-print "Evaluating domain with %d cells" %N
-cell_len = L/N # Origin = 0.0
-points = zeros(N+1,Float)
-
-# Define the centroid points
-for j in range(N+1):
-    points[j] = j*cell_len
+N = 200
 
 # Create domain with centroid points as defined above        
-domain = Domain(points)
+domain = Domain(*interval_mesh(N))
+
 
 # Set initial values of quantities - default to zero
@@ -63 +52 @@
 
 # Set boundry type, order, timestepping method and limiter
-domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4])
+domain.set_boundary({'left': Bd , 'right' : Bd })
+
+
 domain.order = 2
 domain.set_timestepping_method('rk2')
@@ -76 +68 @@
     domain.write_time()
 
-N = float(N)
-HeightC = domain.quantities['height'].centroid_values
+
+exit()
+
+HeightC    = domain.quantities['height'].centroid_values
 DischargeC = domain.quantities['discharge'].centroid_values
-C = domain.centroids
+C          = domain.centroids
 print 'That took %.2f seconds' %(time.time()-t0)
-X = domain.vertices
-HeightQ = domain.quantities['height'].vertex_values
-VelocityQ = domain.quantities['velocity'].vertex_values
-x = X.flat
-z = domain.quantities['elevation'].vertex_values.flat
-stage=HeightQ.flat+z
+X          = domain.vertices
+HeightQ    = domain.quantities['height'].vertex_values
+VelocityQ  = domain.quantities['velocity'].vertex_values
+Z = domain.quantities['elevation'].vertex_values
+Stage      = HeightQ+Z
 
 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
@@ -94 +87 @@
 plot1 = subplot(211)
 
-plot(x,z,x,stage)
+plot(X.flat,Z.flat,X.flat,Stage.flat)
  
 plot1.set_ylim([-1,11])
@@ -101 +94 @@
 legend(('Analytical Solution', 'Numerical Solution'),
            'upper right', shadow=True)
+
 plot2 = subplot(212)
-plot(x,VelocityQ.flat)
+plot(X.flat,VelocityQ.flat)
 plot2.set_ylim([-10,10])
     

anuga_work/development/2010-projects/anuga_1d/channel/test_padarn.py

@@ -2 +2 @@
 import random
 from math import sqrt, pow, pi ,sin, cos
-from channel_domain import *
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
-from config import g, epsilon
+import numpy
+
+from anuga_1d.channel.channel_domain import *
+
+from anuga_1d.config import g, epsilon
 
 
@@ -17 +19 @@
     omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
     t=0.0
-    y = zeros(len(x),Float)
+    y = numpy.zeros(len(x),numpy.float)
     for i in range(len(x)):
         #y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
@@ -26 +28 @@
     N = len(x)
     z_infty = 10.0
-    z = zeros(N,Float)
+    z = numpy.zeros(N,numpy.float)
     L_x = 2500.0
     A0 = 0.5*L_x
@@ -35 +37 @@
 
 def initial_area(x):
-    y = zeros(len(x),Float)
+    y = numpy.zeros(len(x),numpy.float)
     for i in range(len(x)):
         y[i]=initial_stage([x[i]])-bed([x[i]])
@@ -61 +63 @@
 print "Evaluating domain with %d cells" %N
 cell_len = 4*L/N # Origin = 0.0
-points = zeros(N+1,Float)
+points = numpy.zeros(N+1,numpy.float)
 for i in range(N+1):
         points[i] = -2*L +i*cell_len
@@ -95 +97 @@
 
 N = float(N)
-HeightC = domain.quantities['height'].centroid_values
+HeightC    = domain.quantities['height'].centroid_values
 DischargeC = domain.quantities['discharge'].centroid_values
-C = domain.centroids
+C          = domain.centroids
 print 'That took %.2f seconds' %(time.time()-t0)
-X = domain.vertices
-HeightQ = domain.quantities['area'].vertex_values
-VelocityQ = domain.quantities['velocity'].vertex_values
-x = X.flat
-z = domain.quantities['elevation'].vertex_values.flat
-stage=HeightQ.flat+z
+X          = domain.vertices
+HeightQ    = domain.quantities['area'].vertex_values
+VelocityQ  = domain.quantities['velocity'].vertex_values
+Z          = domain.quantities['elevation'].vertex_values
+Stage      = HeightQ + Z
 
 
@@ -113 +114 @@
    
 plot1 = subplot(211)
-
-plot(x,z,x,stage)
+plot(X.flat,Z.flat, X.flat,Stage.flat)
  
 plot1.set_ylim([-1,35])
@@ -122 +122 @@
            'upper right', shadow=True)
 plot2 = subplot(212)
-plot(x,VelocityQ.flat)
+plot(X.flat,VelocityQ.flat)
 plot2.set_ylim([-10,10])
     

anuga_work/development/2010-projects/anuga_1d/generic/generic_domain.py

@@ -104 +104 @@
         self.build_surrogate_neighbour_structure()          
 
-        #Build boundary dictionary mapping (id, edge) to symbolic tags
         #Build boundary dictionary mapping (id, vertex) to symbolic tags
         self.build_boundary_dictionary(boundary)
@@ -745 +744 @@
                 msg += 'in the supplied dictionary.\n'
                 msg += 'The tags are: %s' %self.get_boundary_tags()
-                raise msg
+                raise Exception, msg
 
 
@@ -1247 +1246 @@
                 #print 'name %s j = %f \n'%(name,j)
                 Q = self.quantities[name]
+
                 Q.boundary_values[i] = q[j]
                 #print 'Q=',Q