Changeset 851

Timestamp:

Feb 8, 2005, 7:07:56 PM (20 years ago)

Author:

ole

Message:

First cut at a spatio-temporal boundary.
Need to improve the test, though.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (1 diff)
• generic_boundary_conditions.py (modified) (3 diffs)
• shallow_water.py (modified) (1 diff)
• test_shallow_water.py (modified) (4 diffs)
• util.py (modified) (7 diffs)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -270 +270 @@
                 raise msg
 
+                
     def set_region(self, functions):
         # The order of functions in the list is used.

inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

@@ -187 +187 @@
         #any tagged boundary later on.
         
-        midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
+        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 = {} 
+        
         for i, (vol_id, edge_id) in enumerate(boundary_keys):
             
@@ -199 +203 @@
             if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
             if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])            
-            midpoint_coordinates[i,:] = m
-        
-        self.midpoint_coordinates = midpoint_coordinates    
+            self.midpoint_coordinates[i,:] = m
+            
+            self.boundary_indices[(vol_id, edge_id)] = i
+        
         
-        self.F = file_function(filename, domain, interpolation_points=midpoint_coordinates)
+        self.F = file_function(filename, domain, 
+                               interpolation_points=self.midpoint_coordinates)
         self.domain = domain
 
         #Test
-        q = self.F(0, 0, 0)
+        q = self.F(0, point_id=0)
 
         d = len(domain.conserved_quantities)
-        msg = 'Values specified in file %s must be a list or an array of length %d' %(filename, 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
 
@@ -226 +233 @@
 
         if vol_id is not None and edge_id is not None:
-            x, y = self.midpoint_coordinates[ vol_id, edge_id ]
-            return self.F(t, x, y)
+            i = self.boundary_indices[ vol_id, edge_id ]
+            return self.F(t, point_id = i)
         else:
-            return self.F(t)
+            return self.F(t) #What should the semantics be?
 
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -941 +941 @@
         except Exception, e:
             msg = 'Function %s could not be executed:\n%s' %(f, e)
+            #FIXME: Reconsider this semantics
             raise msg
                     

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -2415 +2415 @@
         Bd = Dirichlet_boundary([0.3,0,0])      
         Bt = Transmissive_boundary(domain1)
-        domain1.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
+        domain1.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
 
         #Initial condition
@@ -2452 +2452 @@
         domain2.reduction = domain1.reduction
         domain2.smooth = False
-        domain2.visualise = True
+        domain2.visualise = False
         domain2.default_order = 2
         
@@ -2465 +2465 @@
         
 
-        #FIXME: Continue from here!!!!!!!
-        return  
-                
         Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
                                       domain2)
-        domain2.set_boundary({'right':Bt, 'bottom':Br, 'diagonal':Bf})
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
 
         #Initial condition 
@@ -2486 +2483 @@
         cv2 = domain2.quantities['stage'].centroid_values
         
-        print take(cv1, (12,14,16))  #Right
-        print take(cv2, (4,6,8))
-        
+        #print take(cv1, (12,14,16))  #Right
+        #print take(cv2, (4,6,8))
         #print take(cv1, (0,6,12))  #Bottom     
-        
-        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
-        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom       
+        #print take(cv2, (0,1,4))
+        #print take(cv1, (0,8,16))  #Diag       
+        #print take(cv2, (0,3,8))       
+        #assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
+        #assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom      
         
         #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)
-    
+        try:
+            os.remove(domain1.filename + '.' + domain1.format)
+        except:
+            pass
+             
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(TestCase,'test')
-    suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
+    suite = unittest.makeSuite(TestCase,'test')
+    #suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
     runner = unittest.TextTestRunner()
     runner.run(suite)

inundation/ga/storm_surge/pyvolution/util.py

@@ -123 +123 @@
     quantities = ['stage', 'xmomentum', 'ymomentum']
     
-    interpolation_points decides at which points interpolated quantities are to be computed
+    interpolation_points decides at which points interpolated 
+    quantities are to be computed whenever object is called.
+    
+    
+    
     If None, return average value
     """
@@ -133 +137 @@
         will be checked and possibly modified
         
-
-
         All times are assumed to be in UTC
         """
@@ -213 +215 @@
         msg = 'File %s must list time as a monotonuosly ' %self.filename
         msg += 'increasing sequence'
-        assert alltrue( time[1:] - time[:-1] > 0 ), msg 
+        assert alltrue(time[1:] - time[:-1] > 0 ), msg  
         
 
@@ -227 +229 @@
             
             ####self.Q = zeros( (len(time),
-             
-            for i, t in enumerate(time[:]):
+            self.T = time[:] 
+            self.index = 0 #Initial index       
+
+            self.values = {}                
+            for name in self.quantities:
+                self.values[name] = zeros( (len(self.T), 
+                                            len(interpolation_points)), 
+                                            Float)                  
+            
+            for i, t in enumerate(self.T):
     
                 #Interpolate quantities at this timestep 
@@ -238 +248 @@
                                          verbose = False)
 
-                for name in self.quantities:                            
-                    print t, name, interpol.interpolate(fid.variables[name][i,:])                    
+                for name in self.quantities:
+                    self.values[name][i, :] =\
+                    interpol.interpolate(fid.variables[name][i,:])   
 
         else:
             pass 
 
-        self.T = time[:]
-        
-        return
-        
-        fid = open(self.filename)        
-        lines = fid.readlines()
-        fid.close()
-        
-        N = len(lines)
-        d = self.number_of_values
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, d), Float)  #Values
-        
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-
-            T[i] = realtime - self.starttime
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-
-
-        self.T = T     #Time
-        self.Q = Q     #Values
-        self.index = 0 #Initial index
-
-        
+    def __repr__(self):
+        return 'File function'
+
+    def __call__(self, t, x=None, y=None, point_id = None):
+        """Evaluate f(t, point_id)
+        
+        Inputs:
+          t: time - Model time (tau = domain.starttime-self.starttime+t)  
+                    must lie within existing timesteps   
+          point_id: index of one of the preprocessed points. 
+                    If point_id is None all preprocessed points are computed
+
+          FIXME: point_id could also be a slice     
+          FIXME: One could allow arbitrary x, y coordinates 
+          (requires computation of a new interpolator) 
+          FIXME: What if x and y are vectors? 
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import zeros, Float
+
+        #Find time tau such that
+        #
+        # domain.starttime + t == self.starttime + tau
+        
+        if self.domain is not None:
+            tau = self.domain.starttime-self.starttime+t
+        else:
+            tau = t
+                
+        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
+              %(self.filename, self.T[0], self.T[1], tau)
+        if tau < self.T[0]: raise msg
+        if tau > self.T[-1]: raise msg        
+
+        oldindex = self.index #Time index
+
+        #Find time slot
+        while tau > self.T[self.index]: self.index += 1
+        while tau < self.T[self.index]: self.index -= 1
+
+        #t is now between index and index+1
+        ratio = (tau - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+                
+        #Compute interpolated values
+        q = zeros( len(self.quantities), Float)
+        
+        for i, name in enumerate(self.quantities):
+            Q = self.values[name]    
+
+            q[i] = Q[self.index, point_id] +\
+                 ratio*(Q[self.index+1, point_id] - Q[self.index, point_id]) 
+
+            
+        #Return vector of interpolated values
+        return q
+                
              
 class File_function_ASCII:
@@ -430 +470 @@
         return 'File function'
 
-        
-
-    def __call__(self, t, x=None, y=None):
+    def __call__(self, t, x=None, y=None, point_id=None):
         """Evaluate f(t,x,y)
 
@@ -438 +476 @@
         result is a vector of same length as x and y
         FIXME: Naaaa
+        
+        FIXME: Rethink semantics when x,y are vectors.
         """