Changeset 623 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Timestamp:

Nov 24, 2004, 4:43:40 PM (20 years ago)

Author:

ole

Message:

Created and tested general File_function (reading and interpolating time series from file) and used it to refactor and simplify File_boundary

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (10 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -518 +518 @@
     """
 
-    #FIXME: Experimental
+    #FIXME: Experimental (from old version). Not in use at the moment
     
     #Shortcuts
@@ -560 +560 @@
             xmomc[k] = ymomc[k] = 0.0
         
+        #FIXME: Delete    
         #From 'newstyle
         #if hc[k] < domain.minimum_allowed_height:        
@@ -901 +902 @@
 def check_forcefield(f):
     """Check that f is either
-    1: a callable object of x,y,t, where x and y are vectors
+    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
@@ -915 +916 @@
         y = ones(3, Float)            
         try:
-            q = f(x, y, 1.0)
+            q = f(1.0, x, y)
         except Exception, e:
             msg = 'Function %s could not be executed:\n%s' %(f, e)
@@ -943 +944 @@
 
 class Wind_stress:
-    """Apply wind stress to water momentum
+    """Apply wind stress to water momentum in terms of
+    wind speed [m/s] and wind direction [degrees]    
     """
 
@@ -959 +961 @@
         map from True north to grid north.
 
-        Arguments can also be Python functions of x,y,t as in
-
-        def speed(x,y,t):
+        Arguments can also be Python functions of t,x,y as in
+
+        def speed(t,x,y):
             ...
             return s
         
-        def angle(x,y,t):
+        def angle(t,x,y):
             ...
             return phi
@@ -978 +980 @@
         forcing_terms as in
 
+        Alternatively, one vector valued function for (speed, angle)
+        can be applied, providing both quantities simultaneously.
+        FIXME: Not yet implemented
+
         domain.forcing_terms.append(W)
-        
         """
 
@@ -1006 +1011 @@
         if callable(self.speed):
             xc = domain.get_centroid_coordinates()            
-            s_vec = self.speed(xc[:,0], xc[:,1], t)
+            s_vec = self.speed(t, xc[:,0], xc[:,1])
         else:
             #Assume s is a scalar
@@ -1019 +1024 @@
         if callable(self.phi):
             xc = domain.get_centroid_coordinates()            
-            phi_vec = self.phi(xc[:,0], xc[:,1], t)
+            phi_vec = self.phi(t, xc[:,0], xc[:,1])
         else:
             #Assume phi is a scalar
@@ -1057 +1062 @@
         ymom_update[k] += S*v        
             
-
-class Wind_stress_from_file:
-    """Apply wind stress read from a file to water momentum
-
-    At this stage the wind field is assumed to depend on time only, i.e
-    no spatial dependency.
-    
-    FIXME: This class may be incorporated in the generic Wind_stress class
-    Superseded
-    """
-
-    
-    def __init__(self, filename):
-        """Initialise windfield from a file with the following format
-
-        time [DD/MM/YY hh:mm:ss], speed [m/s] direction [degrees]
-        
-        e.g.
-        
-        03/09/04 19:15:00, 9.53 39
-        
-        domain.forcing_terms.append(W)
-        
-        """
-
-        import time
-        from Numeric import array
-        from config import time_format
-
-
-        from config import rho_a, rho_w, eta_w
-        self.const = eta_w*rho_a/rho_w
-        
-
-        try:
-            fid = open(filename)
-        except Exception, e:
-            msg = 'Wind stress file %s could not be opened: %s\n'\
-                  %(filename, e)
-            raise msg
-
-
-        line = fid.readline()
-        fid.close()
-        fields = line.split(',')
-        msg = 'File %s must have the format date, values'
-        assert len(fields) == 2, msg
-
-        try:
-            starttime = time.mktime(time.strptime(fields[0], time_format))
-        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] 
-            raise msg
-
-        #Split values
-        values = []
-        for value in fields[1].split():
-            values.append(float(value))
-
-        q = array(values)
-        
-        msg = 'ERROR: File function must return a 1d list or array '
-        assert len(q.shape) == 1, msg
-            
-
-        msg = 'Values specified in file must convert to an array of length 2'
-        assert len(q) == 2, msg
-
-        self.filename = filename
-        self.domain = domain
-        self.starttime = starttime
-
-
-        if domain.starttime is None:
-            domain.starttime = starttime
-        else:
-            msg = 'Start time as specified in domain (%s) is earlier '
-            'than the starttime of file %s: %s'\
-                  %(domain.starttime, self.filename, self.starttime) 
-            if self.starttime > domain.starttime:
-                raise msg
-
-
-        self.read_times() #Now read all times in.
-
-        
-    def read_times(self):
-        from Numeric import zeros, Float, alltrue
-        from config import time_format
-        import time
-        
-        fid = open(self.filename)        
-        lines = fid.readlines()
-        fid.close()
-        
-        N = len(lines)
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, 2), Float)  #Wind field (s, phi)
-        
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            real_time = time.mktime(time.strptime(fields[0], time_format))
-
-            T[i] = real_time - self.start_time
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-        msg = 'File %s must list time as a monotonuosly ' %self.filename
-        msg += 'increasing sequence'
-        assert alltrue( T[1:] - T[:-1] > 0 ), msg
-
-        self.T = T     #Time
-        self.Q = Q     #Windfied
-        self.index = 0 #Initial index
-
-        
-    def __repr__(self):
-        return 'Wind field from file'
-
-        
-
-    def __call__(self, domain):
-        """Evaluate windfield based on values found in domain
-        """
-
-        from math import pi, cos, sin, sqrt
-        
-        xmom_update = domain.quantities['xmomentum'].explicit_update
-        ymom_update = domain.quantities['ymomentum'].explicit_update    
-        
-        N = domain.number_of_elements
-        t = self.domain.time
-        
-        msg = 'Time given in File %s does not match model time'\
-              %self.filename
-        
-        if t < self.T[0]: raise msg
-        if t > self.T[-1]: raise msg        
-
-        oldindex = self.index
-
-        #Find slot
-        while t > self.T[self.index]: self.index += 1
-        while t < self.T[self.index]: self.index -= 1
-
-        #t is now between index and index+1
-        ratio = (t - self.T[self.index])/\
-                (self.T[self.index+1] - self.T[self.index])
-
-        #Compute interpolated values for s and phi
-        q = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
-
-        s = q[0]
-        
-        #Convert to radians
-        phi = q[1]*pi/180
-        
-        #Compute velocity vector (u, v)
-        u = s*cos(phi)
-        v = s*sin(phi)
-
-        #Compute wind stress for this time step
-        S = self.const * sqrt(u**2 + v**2)
-        xmom_update += S*u
-        ymom_update += S*v        
-        
-
-