Changeset 1697 for inundation/ga/storm_surge/pyvolution

Timestamp:

Aug 9, 2005, 5:44:21 PM (20 years ago)

Author:

steve

Message:

Found problem with File_Boundary as used in validation test LWRU2. Have setup new BC Transmissive_Momentum_Set _Stage

Location:

inundation/ga/storm_surge/pyvolution

Files:

• config.py (modified) (1 diff)
• domain.py (modified) (2 diffs)
• generic_boundary_conditions.py (modified) (31 diffs)
• netherlands.py (modified) (3 diffs)
• quantity.py (modified) (10 diffs)
• realtime_visualisation_new.py (modified) (1 diff)
• shallow_water.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/config.py

@@ -87 +87 @@
 if sys.platform == 'win32':
     #default_datadir = 'C:\grohm_output'
-    default_datadir = '.' 
+    default_datadir = '.'
 else:
     #default_datadir = os.path.expanduser('~'+os.sep+'grohm_output')
-    default_datadir = '.' 
+    default_datadir = '.'
 
 

inundation/ga/storm_surge/pyvolution/domain.py

@@ -545 +545 @@
         from config import min_timestep
 
-        timestep = self.timestep
+        # self.timestep is calculated from speed of characteristics
+        # Apply CFL condition here
+        timestep = self.CFL*self.timestep
 
         #Record maximal and minimal values of timestep for reporting
@@ -635 +637 @@
     def update_ghosts(self):
         pass
-
 
     def distribute_to_vertices_and_edges(self):

inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

@@ -10 +10 @@
        A boundary object has one neighbour; the one it
        serves.
-    """   
+    """
 
     def __init__(self):
         pass
-        
+
     def evaluate(self, vol_id=None, edge_id=None):
         msg = 'Generic class Boundary must be subclassed'
@@ -26 +26 @@
     Underlying domain must be specified when boundary is instantiated
     """
-    
+
     def __init__(self, domain = None):
         Boundary.__init__(self)
@@ -33 +33 @@
             msg = 'Domain must be specified for transmissive boundary'
             raise msg
-        
-        self.domain = domain
-        
+
+        self.domain = domain
+
     def __repr__(self):
         return 'Transmissive_boundary(%s)' %self.domain
@@ -43 +43 @@
         of the volume they serve.
         """
-        
+
         q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
-        return q        
+        return q
 
 
@@ -66 +66 @@
     def __repr__(self):
         return 'Dirichlet boundary (%s)' %self.conserved_quantities
-    
+
     def evaluate(self, vol_id=None, edge_id=None):
         return self.conserved_quantities
@@ -88 +88 @@
             msg = 'Function for time boundary could not be executed:\n%s' %e
             raise msg
-        
-
-        from Numeric import array, Float        
+
+
+        from Numeric import array, Float
         try:
             q = array(q).astype(Float)
@@ -96 +96 @@
             msg = 'Return value from time boundary function could '
             msg += 'not be converted into a Numeric array of floats.\n'
-            msg += 'Specified function should return either list or array.' 
+            msg += 'Specified function should return either list or array.'
             raise msg
 
         msg = 'ERROR: Time boundary function must return a 1d list or array '
         assert len(q.shape) == 1, msg
-            
+
         d = len(domain.conserved_quantities)
         msg = 'Return value for f must be a list or an array of length %d' %d
@@ -132 +132 @@
         from config import time_format
         from util import File_function
-        
+
         Boundary.__init__(self)
 
@@ -145 +145 @@
         assert len(q) == d, msg
 
-        
+
     def __repr__(self):
         return 'File boundary'
@@ -174 +174 @@
     File boundary must read and interpolate from *smoothed* version
     as stored in sww and cannot work with the discontinuos triangles.
-    
+
     """
 
@@ -182 +182 @@
         from config import time_format
         from util import file_function
-        
+
         Boundary.__init__(self)
 
         #Get x,y vertex coordinates for all triangles
         V = domain.vertex_coordinates
-        
+
         #Compute midpoint coordinates for all boundary elements
-        #Only a subset may be invoked when boundary is evaluated but
-        #we don't know which ones at this stage since this object can be attached to 
-        #any tagged boundary later on.
+        #Only a subset may be invoked when boundary is evaluated but
+        #we don't know which ones at this stage since this object can be attached to
+        #any tagged boundary later on.
 
         if verbose: print 'Find midpoint coordinates of entire boundary'
         self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
-        boundary_keys = domain.boundary.keys()
-        
-        #Record ordering #FIXME: should this happen in domain.py
-        self.boundary_indices = {} 
-        
+        boundary_keys = domain.boundary.keys()
+
+        #Record ordering #FIXME: should this happen in domain.py
+        self.boundary_indices = {}
+
         for i, (vol_id, edge_id) in enumerate(boundary_keys):
-            
+
             x0 = V[vol_id, 0]; y0 = V[vol_id, 1]
             x1 = V[vol_id, 2]; y1 = V[vol_id, 3]
-            x2 = V[vol_id, 4]; y2 = V[vol_id, 5]            
-
-            #Midpoints        
+            x2 = V[vol_id, 4]; y2 = V[vol_id, 5]
+
+            #Midpoints
             if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
             if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
-            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])            
+            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
             self.midpoint_coordinates[i,:] = m
-            
-            self.boundary_indices[(vol_id, edge_id)] = i
-        
-
-        if verbose: print 'Initialise file_function'
-        self.F = file_function(filename, domain, 
+
+            self.boundary_indices[(vol_id, edge_id)] = i
+
+
+        if verbose: print 'Initialise file_function'
+        self.F = file_function(filename, domain,
                                interpolation_points=self.midpoint_coordinates,
                                verbose = verbose)
@@ -225 +225 @@
 
         d = len(domain.conserved_quantities)
-        msg = 'Values specified in file %s must be ' %filename 
-        msg += ' a list or an array of length %d' %d
+        msg = 'Values specified in file %s must be ' %filename
+        msg += ' a list or an array of length %d' %d
         assert len(q) == d, msg
 
-        
+
     def __repr__(self):
         return 'File boundary'
@@ -247 +247 @@
         else:
             #raise 'Boundary call without point_id not implemented'
-            #FIXME: What should the semantics be?        
+            #FIXME: What should the semantics be?
             return self.F(t)
 
@@ -262 +262 @@
     conserved quantities from a volume as defined by a connection table
     mapping between tuples of (volume id, face id) for volumes that
-    have their boundaries connected. 
+    have their boundaries connected.
 
     FIXME: Perhaps include possibility of mapping between
@@ -274 +274 @@
     def __init__(self, table):
         from domain import Volume
-        
+
         Boundary.__init__(self)
 
@@ -281 +281 @@
 
     def __repr__(self):
-        return 'Connective boundary'        
+        return 'Connective boundary'
 
     #FIXME: IF we ever need to get field_values from connected volume,
@@ -287 +287 @@
     #get_conserved_quantities
     #def get_field_values()
-    
+
     def get_conserved_quantities(self, volume, face=0):
 
         id = volume.id
         if self.connection_table.has_key((id, face)):
-            other_id, other_face = self.connection_table[(id, face)] 
-            
+            other_id, other_face = self.connection_table[(id, face)]
+
             other_volume = self.Volume.instances[other_id]
             cmd = 'q = other_volume.conserved_quantities_face%d' %face;
-            exec(cmd)        
+            exec(cmd)
             return q
         else:
@@ -302 +302 @@
                   %(id, face)
             raise msg
-            
-        
-
-
-
-#FIXME: Add a boundary with a general function of x,y, and t 
+
+
+
+
+
+#FIXME: Add a boundary with a general function of x,y, and t
 
 #FIXME: Add periodic boundaries e.g.:
@@ -317 +317 @@
 #
 #    print i,k
-#    
+#
 #    domain[i].faces[2].neighbour = domain[k].faces[1]
-#    domain[k].faces[1].neighbour = domain[i].faces[2]    
-    
+#    domain[k].faces[1].neighbour = domain[i].faces[2]
+
 
 
@@ -329 +329 @@
     Must specify initial conserved quantities,
     neighbour,
-    domain to get access to model time 
+    domain to get access to model time
     a function f(q_old, neighbours_q, t) which must return
     new conserved quantities q as a function time
@@ -366 +366 @@
         from Numeric import array
         from config import time_format
-        
+
         Boundary.__init__(self)
 
@@ -385 +385 @@
         try:
             starttime = time.mktime(time.strptime(fields[0], time_format))
-        except ValueError:    
+        except ValueError:
             msg = 'First field in file %s must be' %filename
             msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0] 
+            msg += 'I got %s instead.' %fields[0]
             raise msg
 
@@ -397 +397 @@
 
         q = array(values)
-        
+
         msg = 'ERROR: File boundary function must return a 1d list or array '
         assert len(q.shape) == 1, msg
-            
+
         d = len(domain.conserved_quantities)
         msg = 'Values specified in file must be a list or an array of length %d' %d
@@ -414 +414 @@
             msg = 'Start time as specified in domain (%s) is earlier '
             'than the starttime of file %s: %s'\
-                  %(domain.starttime, self.filename, self.starttime) 
+                  %(domain.starttime, self.filename, self.starttime)
             if self.starttime > domain.starttime:
                 raise msg
-        
+
         self.read_time_boundary() #Now read all times in.
-        
-        
+
+
     def read_time_boundary(self):
         from Numeric import zeros, Float, alltrue
         from config import time_format
         import time
-        
-        fid = open(self.filename)        
+
+        fid = open(self.filename)
         lines = fid.readlines()
         fid.close()
-        
+
         N = len(lines)
         d = len(self.domain.conserved_quantities)
@@ -436 +436 @@
         Q = zeros((N, d), Float)
 
-        
+
         for i, line in enumerate(lines):
             fields = line.split(',')
@@ -443 +443 @@
             T[i] = real_time - self.starttime
 
-            
+
             for j, value in enumerate(fields[1].split()):
                 Q[i, j] = float(value)
@@ -449 +449 @@
         msg = 'Time boundary must list time as a monotonuosly '
         msg += 'increasing sequence'
-        
+
         assert alltrue( T[1:] - T[:-1] > 0 ), msg
 
@@ -456 +456 @@
         self.index = 0 #Initial index
 
-        
+
     def __repr__(self):
         return 'File boundary'
@@ -465 +465 @@
 
         t = self.domain.time
-        
+
         msg = 'Time given in File boundary does not match model time'
         if t < self.T[0]: raise msg
-        if t > self.T[-1]: raise msg        
+        if t > self.T[-1]: raise msg
 
         oldindex = self.index
@@ -478 +478 @@
         #if oldindex != self.index:
         #    print 'Changing from %d to %d' %(oldindex, self.index)
-        
-        
+
+
         #t is now between index and index+1
         ratio = (t - self.T[self.index])/\
@@ -485 +485 @@
 
         return self.Q[self.index,:] +\
-               ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
-
-
-    
-    
-        
+               ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
+
+
+
+

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -83 +83 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 100
+N = 50
 
 
@@ -154 +154 @@
 #
 print 'Initial condition'
-domain.set_quantity('stage', Constant_height(Z, 0.0))
-#domain.set_quantity('stage', Constant_stage(-1.0))
+#domain.set_quantity('stage', Constant_height(Z, 0.0))
+domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
 
 #Evolve
@@ -169 +169 @@
 for t in domain.evolve(yieldstep = 0.02, finaltime = 15.0):
     domain.write_time()
+    print domain.quantities['stage'].get_values(location='centroids',indexes=[0])
     #V.update_quantity('stage')
     #rpdb.set_active()

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -182 +182 @@
             self.set_values_from_file(X) #FIXME: More parameters
 
-        
+
         if location == 'vertices' or location == 'unique vertices':
             #Intialise centroid and edge_values
@@ -199 +199 @@
                   Default is "vertices"
 
-        In case of location == 'centroid' the dimension values must
+        In case of location == 'centroids' the dimension values must
         be a list of a Numerical array of length N, N being the number
         of elements. Otherwise it must be of dimension Nx3
@@ -266 +266 @@
         #FIXME: Should check that function returns something sensible and
         #raise a meaningfol exception if it returns None for example
-        
+
         from Numeric import take
 
@@ -404 +404 @@
         #FIXME location, indices
 
-        
+
         import util, least_squares
 
@@ -413 +413 @@
             attribute_name = names[0]
 
-            
+
         if verbose:
             print 'Using attribute %s from file %s' %(attribute_name, filename)
-            print 'Available attributes: %s' %(names) 
-            
+            print 'Available attributes: %s' %(names)
+
 
         try:
@@ -425 +425 @@
                   %(attribute_name, filename)
             raise msg
-        
+
 
         #Fit attributes to mesh
@@ -438 +438 @@
         #Recursively set value using fitted array
         print 'shape', X.shape
-        self.set_values(X, location='vertices')  #FIXME, indexes = None):        
-        
+        self.set_values(X, location='vertices')  #FIXME, indexes = None):
+
 
 
@@ -470 +470 @@
             assert A.shape[0] == len(indexes)
             vertex_list = indexes
-            
+
         #Go through list of unique vertices
         for i_index, unique_vert_id in enumerate(vertex_list):
@@ -606 +606 @@
             #M = 3*m        #Total number of unique vertices
             #V = reshape(array(range(M)).astype(Int), (m,3))
-            
+
             V = self.domain.get_triangles(obj=True)
             #FIXME use get_vertices, when ready
@@ -861 +861 @@
 
             #Gradient
-            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)          
+            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
         else:
             #No true neighbours -

inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

@@ -206 +206 @@
         #print 'update bed image'
         if qname=='elevation':
-            self.vpython_z_models[qname].pos    = self.pos+1.0e-3
-        else:
-            self.vpython_z_models[qname].pos    = self.pos
-
+            self.pos[:,2] = self.pos[:,2]+1.0e-3
+
+        self.vpython_z_models[qname].pos    = self.pos
         self.vpython_z_models[qname].color  = self.colour
         self.vpython_z_models[qname].normal = self.normals

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -113 +113 @@
             from realtime_visualisation_new import Visualiser
             self.visualiser = Visualiser(self,scale_z,rect)
-            self.visualiser.updating['elevation']=True
+            self.visualiser.setup['elevation']=True
         self.visualise = True
 
@@ -181 +181 @@
 
 #         #self.quantities['stage'].interpolate()
-
 
 
@@ -964 +963 @@
 
 
+
+class Transmissive_Momentum_Set_Stage_boundary(Boundary):
+    """Returns same momentum conserved quantities as
+    those present in its neighbour volume. Sets stage
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain = None, function=None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise msg
+
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise msg
+
+        self.domain   = domain
+        self.function = function
+
+    def __repr__(self):
+        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Transmissive Momentum Set Stage boundaries return the edge momentum
+        values of the volume they serve.
+        """
+
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        value = self.function(self.domain.time)
+        q[0] = value[0]
+        return q
+
+
+
+
+
 #class Spatio_temporal_boundary(Boundary):
 #    """The spatio-temporal boundary, reads values for the conserved